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PREFACE 


In  December  1975,  the  Automatic  Data  Processing  (ADP)  Center  of  the  U.  S. 
Army  Engineer  Waterways  Experiment  Station  (WES),  Vicksburg,  Miss.,  submitted  a 
proposal  to  implement  and  evaluate  interval  arithmetic,  a software  system  for 
digital  computer  numerical  analysis,  on  the  Corps  of  Engineers'  primary  engi- 
neering computer — the  WES  Honeywell  G635.  The  proposal  was  later  expanded  to 
include  the  implementation  and  evaluation  of  an  interval  arithmetic  software 
package  on  six  different  computer  systems.  Engineering  and  scientific  data 
problems  were  selected  to  be  used  on  each  of  the  six  computers  with  the  inter- 
val arithmetic  software. 

The  work  was  funded  by  the  Office,  Chief  of  Engineers,  U.  S.  Army,  through 
the  Integrated  Software  Research  and  Development  (ISRAD)  Program,  AT11,  Engi- 
neering Software  Research. 

This  is  Report  3 of  a series  entitled  "Implementation  and  Evaluation  of 
Interval  Arithmetic  Software."  The  other  reports  to  be  published  in  the  series 
are: 

Report  1:  The  State  of  the  Interval:  Evaluation  and  Recommendations 

Report  2:  The  Honeywell  MULTICS  System 

Report  A:  The  IBM  370,  DEC  10,  and  DEC  PDP-11/70  Systems 

Report  5:  The  CDC  CYBER  70  System 

This  report  was  written  by  Mr.  James  Q.  Arnold,  Mr.  Frank  P.  Ford,  and 
Dr.  Richard  G.  Hetherington  of  the  Department  of  Computer  Science,  University 
of  Kansas,  Lawrence,  Kans.  Their  work  was  performed  under  Contract  No,  DACA39- 
76-M-0248,  dated  28  April  1976,  and  Contract  No.  DACA39-77-M-0107,  dated 
24  February  1977,  and  was  partially  supported  by  the  University  of  Kansas 
Computation  Center  project  account  and  by  Research  Grant  3564-5038.  The  work 
concerned  implementation  and  evaluation  of  an  interval  arithmetic  software  sys- 
tem on  the  Honeywell  G635  computer  system. 
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Dr.  J.  Michael  Yohe,  Director  of  Academic  Computer  Services,  University  of 


Wisconsin-Eau  Claire,  developed  and  wrote  the  interval  arithmetic  software 
package  which  was  implemented  on  each  of  the  six  computer  systems.  Dr.  Fred  D. 
Crary,  formerly  with  the  U.  S.  Army  Mathematics  Research  Center,  University  of 
Wisconsin-Madison,  developed  and  wrote  the  AUGMENT  precompiler  which  was  im- 
plemented on  each  computer  system  as  a front-end  to  the  interval  arithemtic 
software  package.  Dr.  Yohe  and  Dr.  Crary  are  specially  thanked  and  recognized 
for  their  technical  contributions  and  assistance. 

Mr.  James  B.  Cheek,  Jr.,  formerly  with  the  ADP  Center,  WES,  provided  ini- 
tial Impetus  and  guidance  for  the  project.  Mr.  Fred  T.  Tracy,  ADP  Center,  WES, 
provided  expert  advice  and  technical  guidance  during  the  project.  Dr.  N. 
Radhakrishnan,  Special  Technical  Assistant,  ADP  Center,  furnished  technical 
guidance  and  general  project  supervision.  The  project  and  the  report  were 
monitored  by  Mr.  William  L.  Boyt  under  the  general  supervision  of  Mr.  D.  L.  Neu- 
mann, Chief  of  the  ADP  Center. 

Directors  of  WES  during  the  project  and  the  preparation  of  the  report  were 
COL  G.  H.  Hilt,  CE,  and  COL  J.  L.  Cannon,  CE.  Technical  Director  was 
Mr.  F.  R.  Brown. 

Copies  of  the  other  reports  of  the  series,  computer  listings  of  the  inter- 
val program  and  of  AUGMENT  for  each  computer  system,  and  runs  of  the  benchmarks 
for  each  computer  system  may  be  obtained  from  the  ADP  Center,  WES. 
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IMPLEMENTATION  OF  THE  'AUGMENT'  PRECOMPILER 
AND  'INTERVAL'  ON  THE  HONEYWELL  G635 

I.  PROJECT  SUMMARY  OF  IMPLEMENTATION  ACTIVITIES 
AUGMENT  version  4K  was  successfully  implemented. 

INTERVAL  II  was  successfully  implemented. 

Test  programs  were  run  against  the  combined  AUGMENT/ INTERVAL  II 
Package . 

Some  Problems  were  identified  and  some  suggestions  are  made. 


II.  INTRODUCTION 


This  report  describes  the  Phase  I activities  of  implementing 
AUGMENT  and  INTERVAL  II  on  the  Honeywell  G635.  As  a result  of  the 
University's  decision  to  release  the  G635  a month  earlier  than 
scheduled,  Phase  I had  to  be  completed  on  the  (compatible)  Honey- 
well 66/60  which  replaced  the  G635.  Unfortunately,  software  errors 
on  both  Honeywell  systems  prevented  completion  of  Phase  I on  schedule. 

Part  III  of  this  report  contains  a discussion  of  the  details 
of  implementing  AUGMENT,  including  a description  of  the  host  software 
errors  mentioned  above,  observations  and  recommendations  regarding 
efficiency,  and  arrangements  for  running  AUGMENT.  Part  IV  contains 
a discussion  of  the  details  of  implementing  INTERVAL  II  in  the  format 
suggested  in  Dr.  Yohe's  draft  report  Implementation  of  the  Package 
on  Other  Hardware,  June  9,  1976.  Part  V contains  a description  of 
the  physical  organization  of  the  AUGMENT/ INTERVAL  II  package  and 
information  related  to  use  of  the  package. 

Listings  of  the  primitives  written  at  K.U.  are  given  in  the 
Appendices,  together  with  results  of  testing  them.  Also,  the 
machine  dependent  constants  needed  by  INTERVAL  II  are  listed  in  an 
Appendix.  Finally,  a number  of  program  listings  together  with  their 
outputs  are  included  in  the  Appendices. 

All  programs  written  at  the  University  of  Kansas  and  all  modi- 
fications to  the  MRC  programming  packages  to  make  them  operational  on 
the  Honeywell  systems  have  been  carefully  and  thoroughly  tested,  and 
are  believed  to  be  correct.  The  University  of  Kansas  and  the  program- 
mers on  this  project  disclaim  responsibility  for  errors  that  may  sub- 
sequently arise  from  use  of  the  programs  and  systems  described  herein. 
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III.  IMPLEMENTATION  OF  AUGMENT 


t 


A.  - :tails  of  bringing  augment  up 

The  initial  effort  to  bring  AUGMENT  up  was  made  on  the  University's 
G635  which  had  200K  words  (36  bits)  for  main  memory,  380  million 
characters  of  disc  storage  and  both  7 and  9 track  tape  drives.  This 
system  supported  batch  and  Time  Sharing  concurrently.  Wherever  pos- 
sible, the  time  Sharing  System  was  employed  in  this  project  in  order 
to  utilize  on-line  editing  capabilities.  For  files  the  size  of  AUGMENT 
this  approach  was  expensive  and  frequently  ran  into  local  boundary 
limits  on  file  space,  processor  time,  and  I/O,... the  typical  large  file 
problems . 

Some  minor  difficulty  was  experienced  in  reading  the  original 
tapes  sent  from  MRC.  In  spite  of  the  obvious  needs  for  standardization, 
it  seems  still  the  case  that  vendor  and  installation  prerogatives 
combine  to  produce  a 'tower  of  Babel'  when  it  comes  to  exchanging  tapes. 
This  was  not  a large  problem,  but  is  noted  here  for  completeness  and 
to  alert  future  users. 

1.  The  Primitives 

After  studying  the  documentation,  it  was  decided  to  rewrite  the 
eight  machine  dependent  primitives  PACK.  CCODE,  M0VH0L,  NUMIN, 
ORDER,  STRCHR,  STRWDS  and  STRLNG.  These  primitives  had  been 
written  for  the  G635  at  the  MRC  and  were  not  in  the  most  efficient 
form.  In  particular,  use  of  such  functions  as  FLD  was  frequent, 
and  produced  grossly  inefficient  object  code  on  the  G635.  All 
eight  primitives  had  been  rewritten  and  were  in  the  final  stages 
of  debugging  (only  one  known  error  still  remained,  albeit  a 
vexing  one)  when  the  G635  was  released  and  the  catastrophic  system 
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errors  on  the  66/60  halted  progress.  Subsequently,  the  versions 
of  these  eight  primitives  which  had  been  sent  from  Wisconsin  were 
re-installed  to  eliminate  all  possible  sources  of  programming 
errors  and  this  is  the  version  currently  being  run.  The  more 
efficient  primitivies  will  be  debugged  as  soon  as  possible  and 
installed  as  a part  of  the  continuing  tuning  of  the  systems.  A 
copy  of  these  routines  will  be  made  available  to  WES  when  they 
have  been  thoroughly  checked  out. 

2.  MIN(J  and  Variable  Dimension  (FDARGS) 

Two  relatively  trivial  errors  appeared  in  the  code  sent  from  MRC. 
MINO  appeared  where  M1N0  was  required  which  was  resolved  by 
having  MINO  call  MIN(5.  Also  non-standard  (for  Honeywell,  at  least) 
use  of  arguments  in  the  routine  FDARGS  produced  a fatal  error. 

Once  identified,  the  correction  was  easy  and  was  checked  with 
Dr.  Crary. 

3.  Compiler  Problems 

During  the  implementation  activity  on  both  Honeywell  systems, 
three  systems  errors  were  identified.  Two  of  the  three  are 
probably  exclusive  with  the  new  version  (SR3I)  on  the  66/60  and 
are  not  errors  on  version  SR8F  for  the  G635.  The  third  error 
resulted  from  a basic  design  flaw  and  exists  on  both  the  G635 
and  the  66/60  (and  all  other  Honeywell  600,  6000,  and  66/ V.'  series 
as  far  as  is  known) . 

In  outputting  error  messages,  a list  of  subroutine  calls 
in  reverse  order  is  provided  by  the  FORTRAN  compiler. 
Intermittantly , strings  of  zeros  are  placed  in  the  name 
field  of  the  list.  In  the  midst  of  other  chaotic  errors, 
this  caused  some  fruitless  searches  for  array  subscript 


systems 
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errors  and  investigations  of  table  construction  algo- 
rithms in  AUGMENT.  Apparently  this  error  does  not  affect 
the  performance  of  AUGMENT,  and  should  be  ignored  if 
it  occurs. 

A much  more  serious  and  disastrously  time-consuming  system 
error  occured  in  the  computation  of  arguments  to  sub- 
programs. After  two  weeks  of  constant  effort  by  the 
entire  project  team,  and  aided  by  frequent  long  telephone 
conversations  with  Dr.  Crary  in  Wisconsin,  it  was  dis- 
covered that  the  66/60  FORTRAN  compiler  was  mishandling 
the  compilation  of  statements  like 

PTR  = SMMAKE (FSTCOL , ENDCOL,  -IABS(CTYP),  K,  0,  -1,  FALSE.) 
(Statements  of  this  form  appear  frequently  in  AUGMENT, 
the  one  given  being  line  AUG  38470  in  PFETCH) . The 
compiler  (version  SR2H  on  the  66/ 60)  produced 

PTR  = SMMAKE (-1ABS (FSTCOL) , K,  0,  -1,. FALSE.) 

Defining  a dummy  variable  for  -IABS(C0L)  and  assigning 

it  prior  to  compiling  the  error  producing  statement  cir- 
cumvented the  problem.  This  error  should  not  appear 
in  SR8F  on  the  G635. 

The  third  error  concerns  the  assembly  language  instruction 
Double  Precision  Floating  Compare  Magnitude  (DFCMG).  If 
DFCMG  is  used  to  compare  two  negative  numbers  whose  mag- 
nitudes differ  by  a large  number  (for  example  more  than 
25 

10  ) the  return  is  always  "greater  than".  In  effect, 

the  argument  in  the  register  is  always  taken  as  having 
magnitude  greater  than  the  argument  in  memory  (even  though, 
for  example,  -1  might  be  in  the  register  and  -10  in 
memory).  This  error  was  first  discovered  on  the  G635  but 
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has  subsequently  been  identified  on  all  large  scale  Honey- 
well systems.  Furthermore,  the  single  precision  form, 
FCMG,  also  falls  as  a result  of  the  same  basic  design 
flaw.  By  Programming,  these  instructions  were  avoided. 
These  errors  are  definitely  part  of  the  G635  system  and 
should  be  publicized  locally. 


4.  Testing  the  Package 

The  test  programs  sent  with  AUGMENT  from  the  MRC  were  run  success- 
fully with  all  modifications  indicated  above. 

One  difficulty  appeared  in  the  way  AUGMENT  handles  the  following 


INTERVAL  A 


PRINT  100,  A(l) , A(2) 

What  happens  is 

(1)  A is  not  taken  over  as  INTERVAL  in  the  AUGMENT  output; 

(2)  the  tables  passed  to  the  host  FORTRAN  are  incorrect 
and  produce  the  fatal  error  message  that  A has  been 
previously  defined  as  a scaler. 

While  it  is  admitted  that  the  PRINT  statement  should  be 
PRINT  100,  A , 

this  is  a potentially  frequent  error.  It  should  produce  an 
AUGMENT  message  and  should  not  place  incorrect  information  in  the 


AUGMENT  output. 


B.  LINKING 

1.  Size  of  the  System 

In  dealing  with  the  space-efficiency  balancing  problem  the  size 
of  the  program  AUGMENT  may  dictate  the  need  for  overlaying.  In 
the  present  implementation  activity,  no  overlays  were  used.  In- 
stead the  emphasis  was  on  efficiency  (especially  run-time  efficiency) 
and  large  amounts  of  file  and  memory  space  were  employed  without 
particular  attention  paid  to  the  benefits  of  overlays.  However, 
in  aa  ongoing-use  environment,  it  is  recommended  that  overlays 
be  employed  to  reduce  costs  and  space  requirements.  This  will  be 
necessary  in  a small  memory  installation;  but  is  desirable  in  any 
case.  An  overlay  version  is  planned  for  implementation  at  K.U. 
and  will  be  available  to  WES. 

2.  Function  Table  Efficiency 

The  current  version  of  AUGMENT  generates  function  tables  for  each 
use.  Considerable  efficiency  might  be  gained  by  establishing 
these  function  tables  once  and  for  all  within  AUGMENT  and  reducing 
the  overhead  costs  to  each  user. 


C.  PLAN  FOR  RUNNING  AUGMENT 


AUGMENT  i3  a precompiler  written  in  FORTRAN  for  FORTRAN 
programs.  It  should  be  maintained  as  an  object  library  on 
tape  or  disk.  It  occupies  approximately  600  blocks  of  disk 
space.  The  library  includes  a mainline  and  two  block  common 
initializing  routines  called  BLDATl  and  BLDAT2.  The  control 
cards  needed  to  access  the  library  are: 


$ 

IDEHT 

$ 

0PTI0N 

FORTRAN,  N0MA1’ 

$ 

F0RTY 

$ 

LIBRARY 

LB 

$ 

EXECUTE 

$ 

FILE 

LB 

$ 

SYS0UT 

21 

$ 

FILE 

20 

*DESCRIBE 

(description  deck) 

*BEGIN 

(FORTRAN  program) 

*END 

$ ENDJOB 

This  setup  reads  a Fortran  program  and  a description  deck 
from  file  05  (the  usual  input  file)  and  outputs  these  on  file  6. 
File  21  is  a copy  of  the  output  of  AUGMENT  which  includes  errors. 
File  20  is  an  S*  file  which  may  be  used  as  an  input  to  the  normal 
Fortran  compiler. 
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IV.  IMPLEMENTATION  OF  INTERVAL 


Due  to  the  later  arrival  from  MRC  of  the  INTERVAL  tape,  and  at 
the  suggestion  of  Dr.  Yohe,  attention  was  first  given  to  the  coding 
of  the  machine  dependent  arithmetic  and  bounding  routines.  After 
reading  the  documents  on  INTERVAL  it  was  thought  that  the  extended 
precision  capability  of  the  G635  could  be  exploited.  Since  all 
floating  point  arithmetic  is  done  in  the  EAQ  registers,  one  has 
available  a full  72  bit  mantissa  which  would  allow  guard  digit 
computations.  However,  a little  experience  with  the  Honeywell 
rounding  options  was  sufficient  to  demonstrate  the  pitfalls  in 
that  direction.  For  one  thing,  the  floating  store  rounded  instruc- 
tion (FST.r)  is  incorrect  in  that  it  is  possible  to  generate  a 
number  and  its  negative,  which,  after  FSTR,  do  not  sum  to  zero. 

FSTR  always  performs  an  "upward"  round.  A second  problem  results 
from  the  fact  that  the  exponent  has  the  same  size  for  both  single 
and  double  precision  which  precludes  use  of  machine  instructions  in 
certain  cases  treated  by  INTERVAL.  As  a result,  the  machine 
rounding  instructions  were  abandoned  and  the  arithmetic  routine 
BPAADD,  BPASUB,  BPAMUL  and  BPADIV  as  well  as  BPACEB  (to  convert 
EXTENDED  to  BPA)  and  BROUND  (to  handle  directed  roundings)  were 
completely  written  in  assembly  language. 

A.  DATA  REPRESENTATION 

Implicit  assumptions  that  determine  choice  of  data  representation 
are  listed  by  Dr.  Yohe  as 

1.  All  operations  and  mathematical  functions  related  to  type  BPA 
will  be  provided  explicitly  in  the  BPA  package. 
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2.  EXTENDED  is  used  in  evaluating  the  BPA  mathematical 
functions,  which  requires  binding  EXTENDED  to  a higher 
precision  data  type  than  CPA. 

3.  A complete  supporting  package  including  all  mathematical 
functions  is  required  for  type  EXTENDED. 

4.  Every  BPA  number  has  an  exact  representation  in  type 
EXTENDED. 

5.  Every  FORTRAN  integer  has  an  exact  representation  in 
type  EXTENDED. 

6.  Conversion  from  REAL  to  BPA  is  exact. 

7.  Bounds  on  the  accuracy  of  mathematical  functions  are 
available. 

In  the  G635,  the  data  type  BPA  has  been  taken  as  identical 
to  REAL,  and  EXTENDED  has  been  taken  as  identical  to  DOUBLE 
PRECISION.  Although  this  greatly  simplifies  the  data  represen- 
tation problems,  items  1,  3 and  7 above  are  not  readily  met  as 
a result  of  deficiencies  in  the  Honeywell  FORTRAN  subroutine 
library  and  lack  of  usable  documentation  for  what  is  there. 

The  following  functions  were  not  supplied  in  the  vendors 
SR8F  FORTRAN  for  the  G635: 


BPA  (REAL)  EXTENDED  (DOUBLE  PRECISION) 


TAN 

DTAN 

DEXP2** 

DTANH 

CBRT 

DARSIN 

DSINH 

DCBRT 

DARCOS 

DCOSH 

DL0G2* 

* Rename  of  DLOG 
**  Rename  of  DEXP2 

I 
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In  the  SR3I  version  of  the  66/60  software,  a new  and  complete 


set  of  FORTRAN  subroutines  are  contained  in  the  FORTRAN  library. 

A copy  of  those  routines  was  requested  from  Honeywell  for  delivery 
with  this  project.  At  the  time  of  this  writing  it  appears  likely 
that  the  request  will  be  grantad.  This  may  be  handled  at  WES  by 
installing  the  new  library  in  place  of  the  current  one  or  by 
treating  it  as  a special  library  for  use  with  the  AUGMENT/ INTERVAL 
package.  In  any  event,  if  the  complete  FORTRAN  library  is  not 
used,  the  C635  FORTRAN  library  will  not  support  INTERVAL  unless 
locally  developed  subroutines  are  made  available. 

B.  CODING  OF  PRIMITIVES 

This  section  contains  a description  of  the  algorithms  developed 
for  the  BPA  arithmetic  and  for  BROUND,  together  with  explanations 
of  the  decisions  made.  Initially  these  algorithms  were  written  to 
incorporate  special  handling  for  the  asymmetric  numbers  described 
below.  However,  after  having  obtained  a working  version  which 
allowed  the  special  cases  to  be  treated  correctly  as  far  as 
INTERVAL  was  concerned,  it  was  decided  that  these  special  cases 
caused  inefficiencies  and  difficulties  of  understanding  beyond 
their  usefulness.  For  this  reason  a change  was  made  to  the  algo- 
rithms which  eliminates  the  special  cases  as  either  inputs  to,  or 
results  from  the  arithmetic  routines.  The  current  versions  con- 
tain these  modifications  as  described  below.  However,  without 
these  cases,  the  requirement  for  passing  arguments  with  the 
correct  sign  to  BROUND  is  no  longer  justified.  This  in  turn, 
implies  that  some  improvement  in  efficiency  will  be  realized  by 


a return  to  end  sign  correction  as  part  of  GROUND  (this  had  been 
abandoned  when  it  was  observed  that  two  passes  through  GROUND  might 
be  required  for  the  special  cases  mentioned  above)  rewriting  all 
the  primitives  in  order  to  accomplsih  this  increase  in  efficiency 
is  a possible  step  in  tuning  the  system. 

Normalized,  floating  point,  two' s-compl emeu t arithmetic  on  the 
Honeywell  635  dones  not  have  symetric  bounds  on  the  numbers  it  can 
represent.  For  single  precision,  normalized  numbers,  the  following 
ranges  are  given; 


i 


negative 

positive 


-(1+2-26) (2_129)>N>-2129 


-129  -27  127 

2 A"sNS(l-2  )(2  ) 


—128 

Furthermore,  zero  is  represented  as  (0) (2  ). 

Two  main  problems  are  caused  by  this  lack  of  symmetry:  the  absolute 
value  of  the  negative  number  with  the  largest  magnitude  cannot  be  repre- 
sented, and  the  smallest,  normalized,  positive  number  has  no  normalized 
negative  counterpart  (although  its  negative  can  be  given  by  an  un-normalized 
number).  To  preserve  the  accuracy  of  the  arithmetic  done  by  INTERVAL,  we 
must  test  for  and  identify  these  cases  at  various  points  in  the  computation. 

Since  the  negative  number  with  the  largest  magnitude  does  not  have 
a negative,  we  do  not  allow  it  to  be  either  an  argument  to  or  an  answer 
from  any  of  the  arithmetic  routines.  If  that  number  is  passed  to  any  of 
the  arithmetic  routines  as  an  argument,  the  overflow  indicator  is  set, 
the  most  negative  number  allowed  by  the  arithmetic  routines  is  loaded 
into  the  A register  of  Interval,  and  control  is  passed  to  BRGUND  to  set 
the  correct  fault  indicator  (overflow  or  infinity)  depending  on  the 
brounding  option  specified  for  that  particular  operation.  If  the  negative 
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number  with  the  largest  magnitude  is  computed  in  any  of  the  arithmetic 
routines,  then  it  is  treated  like  an  overflow  condition,  which  it  actually 
is  to  Interval  - even  though  it  may  not  be  an  overflow  to  the  computer. 

Although  there  is  a negative  number  with  the  same  magnitude  as 
the  smallest  positive  number,  it  is  not  normalized,  and  therefore 
it  cannot  be  returned  as  an  answer  from  the  arithmetic  routines. 
Consequently,  it  is  not  allowed  as  an  argument  either.  If  the 
smallest,  positive,  normalized  number  comes  in  as  an  argument  to 
the  arithmetic  routines,  the  exponent  underflow  indicator  is  turned 
on,  and  control  passes  to  B ROUND . In  BROUND,  it  is  treated  as  a case 
of  underflow  by  one. 

All  arithmetic  routines  pass  a number  with  the  correct  sign  to 
BROUND.  Corrections  for  the  residue  indicator  also  take  place  in  the 
arithmetic  routines. 

//In  these  arithmetic  routines,  MINNEG  is  the  negative 
number  with  the  largest  magnitude,  which  is  representable 
in  the  computer.  NEGBND  is  the  negative  number  with  the 
largest  number  which  is  allowed  in  Interval.  MINPOS  is 
the  smallest,  normalized  positive  number  which  can  be 
represented.  Additional  comments  are  given  for  numbered 
lines .// 

BPASUB:  //compute  A ♦-  A1  - A2  by  negating  A2  and  adding  // 

U + A2 
if  V - 0 

then  A *■  A1 


return 


else  U ♦ -U 

goto  SUB  IN 

fi 

BPAADD:  //compute  A A1  + A2// 

U ♦ A2 
if  U = 0 

then  A A1 
return 

fi 

SUBIN:  A A1 

if  A = 0 

then  A «-  U 
return 

fi 

1 if  Sul  < 1*1 

then  U-*-*-  A 

fi 

if  A < 0 

then  A ■* — A 
U * -U 
SC  «-  1 

fi 

TEMP  *■  E(U)  - E(A) 

2 if  TEMP  > 0 

3 then  E(A)  +■  E(U) 

4 shift  M(A)  right 

fi 


//go  to  add  routine  to  add  AI,A2// 


//if  U • 0,  then  A1  is  the  answer// 


//if  A = 0,  then  the  answer  is  in  U // 


//put  argument  with  smaller// 
//magnitude  in  A// 

//A  must  be  made  positive// 


//we  compute  the  shift  count// 

//if  TEMP  <_  0,  we  need  no  shift// 

TEMP  bits  //line  up  mantissas// 


.•  . . — - 
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‘ 
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15  then  M(A)  «-  M(A)  - (0...01) 

fi 

goto  BROUND 

ft  //end  of  BPAADD// 

Additional  comments: 

lines 

1 At  the  present  time  (1976)  the  DFCMG  instruction  does  not 

work  on  the  Honeywell  635;  consequently,  we  must  work 
around  it,  using  other  instructions  to  do  the  same  thing.  We 
have  done  this  by  making  both  of  the  arguments  positive, 
comparing  them  (which  is  the  same  as  comparing  their 
magnitudes),  and  then  putting  the  smaller  one  in  A and  the 
bigger  one  in  U with  their  correct  signs.  The  problem  with 
the  DFCMG  was  not  known  by  Honeywell  at  the  time  it  was 
reported  to  them. 

2-4  E(U)  - E(A)  will  be  negative  in  one  special  case: 

M(A)  = -M(U)  = 010...0  = %.  A negative  number  normalizes 
to  have  a 1 in  the  first  bit  of  the  mantissa;  consequently, 
its  exponent  is  one  less  than  the  exponent  of  its  absolute 
value.  In  this  case,  we  do  not  want  to  shift  the  exponents. 
When  the  difference  of  the  exponents  is  zero,  even  when  one 
is  negative  and  one  is  positive,  we  do  not  need  to  shift 
the  mantissas,  because  they  line  up  already.  Otherwise, 
we  must  shift  the  mantissa  of  the  A register  right,  to 
line  it  up  with  the  mantissa  of  the  U register.  The 
number  of  bits  we  must  shift  it  is  the  difference, 

E(U)-E(A) . 
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5-10  After  the  addition  of  the  two  arguments,  we  must  restore  the 

sign  of  the  answer.  If  we  produced  an  overflow  in  the 
addition,  simple  negation  will  produce  the  mantissa  with 
the  correct  sign.  One  possible  answer  for  the  overflow 
is  the  smallest  positive  normalized  number.  If  we  negate 
that,  we  will  get  an  underflow,  and  so  after  the  negation, 
we  must  be  sure  to  turn  off  the  EU  indicator,  and  to  turn 
the  EO  indicator  back  on. 

Furthermore,  if  we  produced  an  underflow  in  the  addition, 
one  possible  answer  is  the  negative  number  with  the  largest 
magnitude  representable  on  the  machine.  If  we  try  to  negate 
that,  we  will  get  the  correct  mantissa,  but  we  will  also 
produce  an  overflow  because  that  number  does  not  have  an 
absolute  value  which  is  machine  representable.  Once 
again,  after  we  do  the  negation,  we  must  be  sure  to  turn 
the  EO  indicator  back  off  (if  it  was  turned  on) , and  we 
must  be  sure  to  reset  the  EU  indicator. 

11-15  At  this  point  we  need  to  make  a residue  correction,  and 
we  have  four  cases  to  deal  with: 

1 . SC  = 0 and  A _>  0 

2 . SC  = 0 and  A < 0 

3.  SC  = 1 and  A ^ 0 

4 . SC  = 1 and  A < 0 

Case  1.  Both  A1  and  A2  were  positive.  Thus  the  bits 
lost  decreased  the  magnitude  of  the  answer. 

If  any  of  the  bits  M(A)  28-63  are  ^ ’ the 
correction  has  already  been  made,  so  only  if 
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M(A>2g_£3  = 0 do  we  set  M(A)^^  1. 

Case  2.  A1  and  A2  had  different  signs,  and  the  positive 

had  the  smaller  ma®"itude.  Thus  the  bits 

lost  increased  the  magnitude  of  the  answer. 

If  M(A>2g  = 0,  then  we  must  set  M(A)^  1 

to  correct  for  the  gain  in  magnitude.  If 

M(A)„„  , ^ 0,  however,  the  correction  has  been 
Z o-o  5 

made . 

Thus,  Case  1 and  Case  2 require  identical  treatment. 

Case  3.  A1  and  A2  again  differed  in  sign,  but  this  time 
the  positive  argument  had  the  larger  magnitude. 
Therefore,  the  bits  lost  increased  the 
magnitude  of  the  answer.  If  28-63  ^ ® 
and  A > 0,  any  truncation  of  non-zero  bits  in 
M(A)_0  ,,  will  insure  proper  brounding  when 
augmentation  is  implied.  If  28-63  = 
however,  truncation  will  leave  the  answer  too 
large,  and  so  we  must  subtract  0...01  from 
M(A)  to  insure  proper  brounding  when  truncation 
is  implied. 

Case  4.  Both  A1  and  A2  were  negative,  and  thus  the 
non- zero  bits  lost  decreased  the  magnitude 
of  the  answer.  If  MCA^g.gj  ^ 0,  then 
truncation  will  insure  correct  augmentation, 
while  the  non-zero  bits  will  guarantee  proper 
rounding  towards  zero.  But  once  again, 
if  M(A) 2g_g 3 = °»  neither  augmentation  nor 
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brounding  toward  zero  will  work,  and  so  we 


i 


1 


j 

i 


must  subtract  0...01  from  M(A) . 

Thus  Cases  3 and  U require  identical  treatment. 

BPAMUL,:  //compute  A-<-Al*A2// 

//convert  arguments  A1,A2  to  positive  values.  Set  SC:=1 
if  answer  must  be  sign  - corrected  before  bround// 
if  A1  < 0 

then  A1+— Al 
SC<-1 

fi 

if  A2<0 

then  A2+— A2 
SC-^l-SC 

fi- 

ZZ  compute  A1*A2  without  normalizing// 

A+A 1 *A2 

//handle  overflow  or  underflow  from  above// 
if  exponent  overflow  then  go  to  BPAME0  fi 
if  exponent  underflow  then  go  to  BPAMEU  fi 
Abnormalized  A 

if  exponent  underflow  from  normalization  then  go  to  BPAMEU  fi 
if  SC=1  then  A-*— A fi 
go  to  BROUND 

BPAMEU : A«--A 

0PTI0N+rounding  option 
BPAFLT  b 0 
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if  SOI  then  A*— A fi 
go  to  .EU. 

// .EU.  is  a section  of  BR0UND  for  exponent  underflow// 


BPAME0: 


BPAD1V : 


I 


0PTI0N  *-  rounding  option 
BPAFLT  +0 

//normalize  A.  This  may  fix  exponent  overflow// 

E(A)=0 

if  SC=1  then  A ♦— A fi 
Abnormalized  A 
//E(A)  has  shift  amount// 

E(A)  = E(A)  + E(A1)  + E(A2) 
if  exponent  overflow  then  go  to  .E0.  fi 
//compute  A-*-Al/A2// 

// convert  to  positive  values// 
if  A1  < 0 then 
Alb-Al 
SOI 

fi 

if  A2=0  then  go  to  ZER0DV  fi 
if  A2  <0  then 
A2-<-_A2 
SOl-SC  fi 

//determine  if  dividend  larger  than  divisor.  This  section 

(specially  DVF  instruction)  works  only  if  dividend  < divisor.// 
if  M(A1)>M(A2)  then 
M(A)  <-  M(A)  /2 
SHIFT- 1 

fi 
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if  M(A1)  = M(A2) 
SHIFT  - 0 


EXP0: 


EXPU : 


fi 

A«-  quotient  (M( A1 ) /M(A2 ) ) 

Q+-  remainder  (M(A1 ) /M(A2) ) 
if  Q + 0 then  RI:»1  fi 
//RI  is  residue  indicator// 

//negate  if  necessary.  Clear  exponent  to  avoid  its  interference// 
E(A)=0 

if  SC  = 1 then  A-*~A  fi 

//now  compute  exponent.  (A)  is  from  normalization// 

E(A)  - SHIFT 
A *■  normalized  A 
E(A)  = E(A)  + E(A1)  - E(A2) 
if  e(A)  > 128  then  go  to  EXP0  fi 

if  E(A)  <-128  then  go  to  EXPl'  fi 
if  RI  = 1 then  M(A)=  M(A)+l  fi 
dPTI(/N  -<-0 
BPAFLT  +-0 

go  to  . E0 . 

0PTI0N  rounding  option 
BPAFLT  <-0 
go  to  . EU . 


Understanding  the  brounding  of  two's  complement  numbers  requires 
some  familiarity  with  the  way  those  numbers  are  represented,  particularly 


the  negative  numbers. 

A positive,  single-precision,  normalized,  floating-point  number 

in  the  Honeywell  635  must  have  a 0 in  bit  zero,  and  a 1 in  bit  one. 

After  that,  it  may  have  any  combination  of  l's  and  0's.  Consequently, 

the  mantissa  may  represent  a number  as  small  as  2 *,  and  as  large  as 
-27 

1-2  . If  we  add  bits  to  a positive  mantissa  or  concatenate  bits 

at  the  right  end  of  a positive  mantissa,  and  those  bits  are  nonzero, 
then  we  increase  the  magnitude  of  the  mantissa;  if  we  drop  a non-zero 
bit  off  the  right  end  of  a positive  mantissa,  we  will  decrease  the 
magnitude  of  th2  number  represented  by  the  mantissa. 

The  Honeywell  635  has  a 28-bit  single-precision  mantissa,  and  a 
64-bit  double-precision  mantissa.  The  rounding  options  for  Interval 
are  accomplished  as  follows: 

U,A:  If  any  of  the  bits  in  M(A) 28-63  are  non-zero»  c^en 
we  want  to  return  the  next  higher  single-precision 
mantissa  as  the  answer.  If  M(A)  28-63  = tllougb, 
we  don't  want  to  change  the  answer  we  have.  Thus 
we  want  to  add  a number  to  A that  will  propagate 
a carry  from  any  of  the  bits  28-63  into  bit  27  of 
the  mantissa  without  generating  a carry  if  M(A)2g  “ 0. 

L,T:  Any  non-zero  bits  in  M(A) 28-53  "^ke  the  answer  too 


large,  so  we  truncate  them  and  use  only  the  single 
precision  answer.  If  the  extra  precision  word  is 
all  zero,  then  we  may  simply  truncate  In  that  case  too. 


up,  while  a 0 means  we  round  down.  Thus,  we  may  aa  a 
one  into  M(A) 2g  and  we  will  get  the  correct  answer. 

If  there  is  a one  there,  we  will  generate  a carry 
into  the  single  precision  word,  if  there  is  a zero, 
then  no  carry  will  be  generated,  and  the  answer 
will  be  rounded  down  as  desired. 


Negative  numbers  are  somewhat  different.  First,  they  normalize 

differently;  they  must  have  a 1 in  the  first  bit,  a zero  in  the  second 

bit,  and  after  that,  any  combination  of  l's  and  0's  may  occur.  Thus, 

-1  -27 

the  range  of  negative  mantissas  is  from  -1  to  -(2  +2  ).  A negative 

mantissa  may  represent  a value  with  a greater  magnitude  than  a positive 
one,  but  a positive  mantissa  may  represent  a value  with  a smaller  magnitude 
than  a negative  one.  If  we  add  bits  to  a negative  mantissa  or  concatenate 
bits  at  the  right  end  of  a negative  mantissa  and  those  bits  are  non-zero, 
we  will  decrease  its  magnitude;  if  we  chop  a non-zero  bit  off  the 
right  end  of  a negative  mantissa,  we  will  increase  the  magnitude  of  the 
number  represented  by  the  mantissa.  This  is  very  important  to  our 
b rounding.  The  brounding  options  are  executed  as  follows: 

U,T:  Both  U and  T imply  brounding  toward  zero;  if  any 

of  the  bits  in  MCA^g^-j  are  one,  then  we  want  to 
generate  a carry  into  M(A)q  ^7  to  decrease  the 
magnitude  of  the  answer.  Only  if  M(A)2g_g^  = ® 
will  we  not  want  to  generate  that  carry.  Thus 
we  may  add  the  same  constant  here  as  we  did  for 
augmentation  in  the  positive  case. 

L,A:  To  find  the  lower  bound,  any  ones  in  MCA^g^^  must 
be  ignored,  because  they  decrease  the  magnitude  of 


the  answer,  and  to  find  the  lower  bound  (or  to 


augment)  we  must  Increase  the  magnitude.  If  there 
are  no  l's  there,  we  have  the  correct  answer.  Thus 
for  both  L and  A options  here,  we  may  simply  chop 
off  the  double  precision  word. 

R:  If  we  are  rounding  a negative  number,  then  we  want 

to  go  away  from  zero  if  the  extra  precision  word 
has  a 1 in  bit  28  and  0's  in  M(A) 29-63  or  if  bit 
28  has  a zero  in  it.  To  take  the  answer  away  from 
zero,  we  want  to  truncate  those  bits.  On  the  other 
hand,  if  M(A)7g  = 1,  and  any  of  M(^^29-63  = f*1611 

we  want  to  bround  toward  zero.  Consequently,  we  want 
to  generate  a carry  into  MiA^-j  of  M(A)^g  = 1 and 
any  of  the  bits  in  MiA^^  gj  = 1 . 

One  important  thing  to  realize  about  this  brounding  is  that 
overflow  and  underflow  may  occur  in  certain  conditions.  We  may  get 
an  overflow  when  augmenting  positive  numbers, 'and  we  may  get  an  under- 
flow by  one  if  we  truncate  negative  numbers.  These  conditions  must 
be  checked  when  the  brounding  is  done,  and  the  appropriate  steps 
must  be  taken  if  a fault  is  recognized. 

Three  numbers  are  used  to  do  brounding:  UTA  is  used  for  the  U 
option  (positive  and  negative) , the  T option  (negative),  and  the  A 
option  (positive);  NEGRND  is  used  by  the  R option  for  negative  numbers; 
POSRND  is  used  by  the  R option  for  positive  numbers.  The  mantissas 
of  the  three  numbers  are  presented  here. 


UTA 


1 00. 

. .00 

1 11.. 

Till 

0 

27 

28 

63 

[ 00.  , 

. .00  1 

[ Oil. 

..11  | 

0 

27 

28 

63 

| 00., 

. .00  ] 

[ 100. 

..00  | 

0 27  28  63 

When  one  of  these  numbers  is  used,  its  exponent  is  set  to  the  exponent 
of  the  number  being  brounded. 

Finally,  we  may  present  the  following  decision  table  for  the 
brounding  routine. 


+ 

- 

u 

add  UTA 

add  UTA 

L 

none 

none 

T 

none 

add  UTA 

A 

add  UTA 

none 

R 

add  POSRNI 

add  NEGRND 

In  the  following  brounding  routine,  several  variables  are  used, 
and  they  represent  certain  floating-point  constants. 

MINNEG  is  the  most  negative  machine  representable  number. 

NEGBND  is  the  most  negative  number  allowed  as  an  argument 
to  or  as  an  answer  from  the  arithmetic  routines. 

MAXNEG  is  the  negative,  single-precision,  normalized 
number  closest  to  zero. 

MINPOS  is  the  smallest,  normalized,  positive  number 
which  can  be  represented  in  the  machine. 

POSBND  is  the  smallest,  single  precision,  normalized 

number  allowed  as  an  argument  or  an  answer  from  the 
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arithmetic  routines 


MAXPOS  is  the  largest,  positive  number. 

UTA,  NEGRND,  and  POSRND  are  as  described. 

B ROUND:  normalize  A 

if  CO  + 'U'  and  CO  t 'L'  and  CO  t 'T'  and  CO  * 'A'  and  CO  ^ 'R' 

then  CO  •+■  'R'  //if  correction  option  is  invalid,  use  R/ / 

fi 

case 
:E0  = 1: 

if  K >_0 

then  A MAXPOS 

if  CO  = 'A'  or  CO  « 'U'  or  CO  = 'R' 
then  INF  1 

fi 

returr 

else  A f-  NEGBVD 

if  CO  = 'A'  or  CO  = 'L'  or  CO  = 'R' 
then  INF  *■  1 

fi 

return 

ft 

:EU  - 1: 

if  A > 0 

then  if  CO  = 'A'  or  CO  = 'U'  or  (CO  = 'R'  and  underflow  by  1) 
then  A POSBND 
else  A «-  0 

fi- 
re turn 
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//  find  the  upper  bound  on  A// 


//overflow  may  occur  when  A is  positive// 


//underflow  may  occur  when  A is  negative// 


else  if  CO  = 'A'  or  CO  = 'L'  or  (CO  = 'R'  and  underflow  by  1) 
then  A * MAXNEG 
else  A *-  0 

fi 

return 

fi 

:C0  = ’U' : 

A * A + UTA 
if  EO  = 1 

then  INF  1 

A +-  MAXPOS 
else  if  F.U  = 1 

then  A *■  0 

fi 
fi- 
re turn 
: CO  = 'L' : 

return 
:C0  = 'T' : 
if  A > 0 

then  return 

else  A +■  A + UTA  //A  < 0// 
if  EU  = 1 

then  A *■  0 

fi 

return 

fi 

:C0  = 'R':  //round  to  nearest  machine  number// 

if  A > 0 


//find  lower  bound// 

//store  single-precision  in  all  cases// 
//b round  toward  zero// 

//store  single-precision  if  positive// 


-w  ■ ..  , 
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then  A * A + POSRND 


if  EO  = 1 

then  A MAXPOS 
INF  * 1 

fi 

return 

else  A <-  A + NEGRND 

£/  EU  = 1 

then  A «-  MAXNEG 

fi 

return 


ft 


: CO  = ’A’: 
vf  A > 0 

then  A «-  A + UTA 
if  EO  = 1 

then  A ■*-  MAXPOS 
INF  i-  1 

fi- 
re turn 
else  return 

fi 


end 

//end  of  bround// 
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//A  < 0// 


//augment  A// 


//A  < 0// 


C.  MACHINE  DEPENDENT  CONSTANTS 


The  machine  dependent  constants,  given  in  BPA  format,  which 
must  be  provided  to  INTERVAL  are  listed  in  fl].  These  constants 
are  located  in  COMMON/ BP ACON/  and  for  the  G635  version  the  actual 
numbers  used  are  given  in  Appendix  B.  iifhile  many  of  these  re- 
quired constants  are  obtained  in  a straightforward  manner 
I 

(e.g.  the  left  endpoint  of  the  INTERVAL  containing  PI),  several 
must  be  derived  from  vendor  documentation  (e.g.  half-length  of 
INTERVAL  about  1 where  LN  accuracy  decreases).  The  approach 
taken  in  determining  these  constants  was  the  pragmatic  one  of 
1 comparing  output  from  the  routines  with  published  tables 

r 2, 3, 4, 5, 6].  It  is  expected  that  some  improvement  in  the  INTERVAL 
results  might  be  realized  by  spending  more  time  on  the  evaluation 
of  these  constants. 

One  additional  comment  in  passing:  the  description  of  the 
constant  FRACBD  does  not  clearly  bring  out  the  fact  that  FRACBD  is 
that  constant  which  when  added  to  any  number  effectively  truncates 


the  decimal  part  and  produces  a real  integer. 


D.  ERROR  BOUNDS  FOR  EXTENDED  PRECISION  ROUTINES 


As  described  in  fl],  the  method  of  implementing  interval 
mathematical  functions  is  to  assign  an  error  term  to  each 
function  which  is  then  utilized  in  the  obvious  way  to  determine 
the  value  (interval)  taken  by  the  function.  Let  f be  a con- 
tinuous interval  valued  function  of  an  interval  variable  which 
we  write  as 

f ( T a . b ] ) = Tc.d]. 

Then,  there  exist  points  a'  and  b'  in  [a,b]  such  that  f(a')  = c 

and  f(b')  = d.  If  the  EXTENDED  library  routine  (double  precision 

routine  in  the  present  instance)  produces  the  value  f_  with  an 

h 

error  bound  €,  then  the  result  computed  by  INTERVAL  II  is 

£([a,b])  = fV(fE  (a')-€,  A (f^  (b’)+€)] 
where  V and  A are  downward  and  upward  directed  roundings, 
respectively. 

From  this  formulation,  it  is  seen  that  accurate  error  bounds 
are  essential  to  optimize  interval  function  evaluations. 

However,  it  is  extremely  difficult  to  calculate  these  bounds 
accurately  and  hence  the  conservative  approach  suggested  by 
Dr.  Yohe  in  Tl]  has  been  adopted.  In  effect  this  means  using 
the  error  bounds  specified  in  the  vendor  documentation  for  these 
routines. 


E.  MODIFICATIONS  TO  INTRAP 


I 


INTRAP  provides  a printed  record  of  faults  encountered 
in  interval  operations.  It  is  assumed  that  each  routine 
which  can  call  INTRAP  has  at  most  three  arguments,  and 
that  the  following  information  is  provided  to  INTRAP: 

ID  a 3-character  suffix  specifying  the  calling  routine 

TYPA  the  type  of  argument  A 

TYPB  the  type  of  argument  B 

TYPR  the  type  of  the  result 

Values  of  TYPA,  TYPB,  TYPR  may  be  any  one  of 

0 null 

1 FORTRAN  INTEGER 

2 FORTRAN  REAL 

3 FORTRAN  DOUBLE  PRECISION 

4 BP  A 

5 EXTENDED 

6 INTERVAL 

Because  BPA  was  identified  with  REAL  and  EXTENDED  was 
identified  with  DOUBLE  PRECISION,  the  built-in  conversion 
routines  were  suitable  for  this  implementation  and  no 
new  code  was  added  to  INTRAP . 

F.  PROCESSING  WITH  AUGMENT 

INTERVAL  consists  of  a description  deck  for  use  with 
AUGMENT  and  a package  of  compatible  subroutines  and  functions 
to  handle  INTERVAL  data.  It  should  be  maintained  as  two 
files — a library  file  and  a description  deck.  The  control 


35 


cards  needed  for  it  would  be: 


$ IDENT 

. (AUGMENT  CONTROL  CARDS) 

$ FILE  20,  AISD 

* DESCRIBE 

$ SELECT  INTERVAL  DECK 

* BEGIN 

Fortran  program 

* END 


$ 

OPTION 

FORTRAN,  NOMAP 

$ 

FORTY 

$ 

LIBRARY 

LC 

$ 

FILE 

S* , AIDD 

$ 

EXECUTE 

$ 

PRMFL 

LC,  INTERVAL  LIB 

Fortran 

data 

$ 

ENDJOB 

»• 


G.  CHECKING  THE  PACKAGE 

The  test  decks  supplied  by  Dr.  Yohe  were  run  against 
INTERVAL  II  in  the  manner  described  in  the  preceding  section, 
and  the  output  was  checked  with  the  sample  output.  As  mentioned 
in  section  E,  these  tests  utilized  INTRAP  as  it  was  sent 
to  us  by  MRC.  The  machine  dependent  constants  described  in 
C and  D were  used  in  these  tests.  The  results  matched 
and  the  test  was  assumed  passed.  Appendix  C contains  these 
test  results. 


H.  TUNING  THE  PACKAGE 


Tuning  the  package  is  interpreted  to  mean  taking 
any  modifications  to  code  or  data  which  will  tesult  in 
reducing  time  or  space  requirements,  or  will  produce  sire 
accurate  interval  results.  (Conflicts  may  occur  between 
these  efforts;  for  example,  run  time  may  be  increased  in 
order  to  get  tighter  intervals.)  The  major  activities 
planned  for  tuning  the  system  are: 

1.  Replacing  AUGMENT  Primivitives  discussed  in  section 
IIIA.l.  of  this  report. 

2.  Constructing  an  overlaved  version  of  the  package. 

(B-l  above) 

3.  Fixing  the  AUGMENT  function  tables  rather  than 
regenerating  them  for  each  use.  (B-2  above) 

4.  Improve  error  bounds  for  extended  precision  routines 
(C  and  D above) 

5.  Remove  all  unnecessary  calls  in  the  AUGMENT  output 
(as  suggested  in  [1]) 

6.  Rewrite  arithmetic  primitives  for  INTERVAL  as  discussed 
in  III.B  of  this  report. 


V. 


USING  THE  COMBINED  AUGMENT/ INTERVAL  PACKAGE 


t 


A.  ORGANIZATION  (PHYSICAL  DESCRIPTION) 

The  particular  organization  to  be  used  at  WES  will  depend  to 
great  extent  on  local  decisions.  The  modules  being  delivered  include 

1.  Object  library  of  AUGMENT  programs. 

2.  Source  library  of  AUGMENT  programs. 

3.  Object  library  of  INTERVAL  programs. 

4.  INTERVAL  description  deck  for  AUGMENT. 

A perm  file  containing  control  cards  (see  III.C)  needed  to 
load  and  execute  the  AUGMENT  object  library  and  a perm  file  con- 
taining control  cards  (see  IV. F)  to  run  AUGMENT  and  the  description 
deck  will  be  created  and  accessed  by  $ SELECT  cards.  Another  perm 
file  to  describe  the  necessary  libraries  and  the  output  file  from 
AUGMENT  to  permit  the  FORTY  compiler  to  execute  the  AUGMENT  output 
will  also  be  created  and  accessed  by  a $ SELECT  card. 

B.  USER  INFORMATION 

For  input/output,  it  should  be  noted  that  a variable  of 
type  INTERVAL  is  converted  to  a variable  of  type  REAL  with  dimen- 
sion 2.  Therefore,  normal  FORTRAN  I/O  methods  may  be  used. 

Section  IV  F contains  specific  control  card  information  for 
running  programs. 


! 
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EVALUATION  OF  INTERVAL  ARITHMETIC 
FOR  THE  HONEYWELL  G635 


VI.  INTRODUCTION 

At  its  instigation  this  project  was  viewed  as  a relatively  routine 
implementation  of  the  Augument/Interval  programming  system  on  the  GE/ 
Honeywell  635  computer  followed  by  a set  of  benchmark  runs  to  investigate 
the  efficacy  of  using  Interval  Arithmetic  to  analyze  error  in  the  problem 
solution  methods  for  use  at  WES.  From  the  beginning,  difficulties  were 
encountered  in  all  aspects  of  the  project  activity,  including  operating 
system  failures  and  misunderstood  communications  at  K.U.;  delays  in 
delivery  of  tapes,  inefficiencies,  and  occassional  errors  in  coding  sent 
by  MRC;  and  untested  benchmark  programs  and  data  sent  by  WES.  While  some 
difficulties  of  this  sort  are  to  be  expected,  the  sum  total  of  all  of 
them  coupled  with  the  rather  poor  performance  of  the  initial  untuned 
Interval/Augument  system  lead  to  extensive  delays  in  the  completion  of 
the  contracted  work.  Every  effort  has  been  made  to  satisfy  the  letter 
and  the  spirit  of  the  contractual  agreement  on  the  parts  of  all  parties 
involved.  However,  Che  simple  fact  is  that  the  task  was  monumentally 
greater  than  anticipated  and  even  now  leaves  many  areas  of  potentially 
fruitful  further  investigation  unexplored. 

Briefly,  the  chronology  of  events  since  the  report  on  implementation 
given  November  1976  at  Vicksburg  is: 

Corrected  versions  of  the  benchmark  algorithms  were 
received 

Overlay  version  of  AUGMENT/ INTERVAL  implemented 

CAUSSE  benchmark  runs  initiated  after  receiving 
the  complete  data 

SPLINE  benchmark  run  successfully 
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Nov,  - Dec.  1976 
January  1977 

February  1977 
March  1977 
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Report  on  completion  of  Phase  III  sent  to  WES 

FFT  benchmark  runs  initiated 

INTERVAL  arithmetic  and  brounding  routines 
rewritten 

BPA  subroutine  calls  eliminated 

Reruns  of  benchmarks  for  timing  comparisons 


March  1977 
April  1977 

May  - June  1977 
July  1977 
July  1977 


■ 


VII.  DISCUSSION  OF  DIFFICULTIES  ENCOUNTERED  IN  USING  INTERVAL  ARITHMETIC 


A.  DEPENDENCY 

There  is  inherent  in  the  operations  of  interval  arithmetic  a defect 

which  may  increase  width  of  intervals  in  computed  results  unrealistically. 

The  generic  term  describing  this  defect  is  dependency.  It  derives  from 

endpoint  calculations  which  make  use  of  expressions  which  are  dependent, 

but  the  interval  arithmetic  operations  treat  them  as  independent.  Thus, 

X 

if  Y = — , where  X = [1/2,2],  the  definition  of  interval  division  produces 
Y = [^-,4].  This  result  includes  the  desired  (correct)  answer  [1,1]  but 
at  the  same  time  has  increased  the  interval  width  substantially.  Other 
anoroolies  include  the  fact  that  additive  and  multiplicative  inverses 
don't  exist  in  the  expected  forms:  if  X = [a,b],  a b,  then  X - X = 

[a  - b,  b - a]  which  is  not  [0,0]  (unless  a = b)  and  X/X  = [r  , ~]  1*  [1,1] 

D 3 

(unless  a = b ^ 0) . The  problem  is  that  the  two  occurrences  of  the  same 

interval  X,  in  these  illustrations,  are  treated  by  the  interval  arith- 

metic operations  as  independent  when  in  fact  they  are  dependent.  There 
is  no  simple  resolution  to  this  difficulty  because  the  definitions  of 
the  interval-arithmetic  operations  must  treat  the  two  interval  operands 
as  independent  in  order  to  guarantee  inclusion  of  the  correct  result  in 
the  interval  answer.  Program  testing  of  the  operands  for  dependency  is 
hopelessly  complicated  in  any  but  the  simplest  sequences  of  operations. 

The  fact  is  that  interval  arithmetic  incorporates  an  error  generation 
characteristic  which  is  distinct  from  ordinary  arithmetic  error  considera- 
tions. The  result  is  that  proven  methods  and  algorithms  of  real 
arithmetic  often  fare  badly  when  converted  directly  into  interval 


arithmetic . 
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The  dependency  problem  is  most  effectively  dealt  with  by  avoidance  - 


simply  not  applying  interval  analysis  to  those  algorithms  or  data  which 
are  known  to  produce  large  dependency  error.  Unfortunately,  there  do  not 
yet  exist  practical  methods  for  deciding,  before  the  fact,  what  will  be 
the  efficacious  interval  algorithms  and  what  will  not.  Experimentation 
often  is  the  most  direct  means  of  obtaining  such  information  although  some 
guidelines  have  been  developed  by  F.N.  Ris  [13]  and  E.R.  Hansen  [12].  The 
upshot  of  this  situation  for  an  installation  such  as  WES  is  that  many 
production  programs  already  in  use  cannot  really  be  effectively  analyzed 
as  interval  arithmetic  programs.  In  those  instances  where  'good'  results 
are  obtained  from  the  interval  runs,  one  can  be  reasonably  sure  that  the 
real  arithmetic  algorithm  is  numerically  stable.  But  when  interval  results 
are  meaningless  because  the  intervals  are  too  wide,  then  one  cannot  conclude 
that  the  real  arithmetic  algorithm  was  also  unstable  - only  that  the  increase 
in  interval  widths  was  too  great.  If  the  dependency  width  (a  term  used  by 
F.N.  Ris  [13]  to  describe  that  part  of  the  increased  interval  width  due 
exclusively  to  dependency)  is  the  principal  contributor  to  the  increased 
widths  of  interval  results,  then  the  algorithm  should  not  be  analyzed 
using  interval  arithmetic  and  a different  analysis  technique  should  be 
employed.  On  the  other  hand,  if  dependency  width  is  small,  even  though  the 
interval  width  of  results  is  too  large,  the  algorithm  may  be  salvaged  by 
some  one  or  more  of  the  traditional  real  arithmetic  error  reduction  pro- 
cedures, e.g.,  using  higher  precision  calculations,  more  terms  in  the 
approximations,  etc.  Of  course,  such  modifications  to  an  algorithm  may 
cause  it  to  exhibit  dependency  which  it  previously  did  not  have.  To  illustrate, 
consider  the  simple  series 


£(x) 


, ^ 2 3 

-1  + X - X + X 


4 

- X ... 

If  the  first  two  terms  are  used,  there  is  no  dependency  problem.  But 
should  two  terms  provide  insufficient  accuracy,  adding  the  third  term 
will  improve  accuracy  and  at  the  same  time  introduce  a dependency  problem. 

B.  COMPARISON  OF  INTERVALS 

Dependency  is  by  far  the  overriding  consideration  in  using  interval 

arithmetic.  If  significant  in  an  algorithm,  it  may  obliterate  any 

meaningful  interval  results.  However,  other  sources  of  difficulty  may 

enter  the  interval  arithmetic  operations  for  a particular  solution 

method.  One  of  these  is  the  problem  encountered  in  translating  the 

compare  operation  from  real  arithmetic  to  interval  arithmetic.  An 

intuitive  approach  might  be  to  compare  endpoints  and  when  the  sup  of 

one  interval  is  less  than  the  inf  of  the  other  interval,  then  say  the 

first  interval  is  less  (smaller)  than  the  second.  This  works  fine  for 

non-overlapping  intervals.  When  the  intervals  overlap,  an  appeal  might 

be  made  to  a real  comparison  of  their  midpoints.  But  saying  that  the 
12  51 

interval  [-^j  , Yqq]  is  less  than  the  interval  [0,1],  as  would  follow, 

leads  to  other  difficulties  since  nearly  half  the  numbers  in  [0,1]  are 

12  51 

less  than  all  values  in  [-jjr  , Yoo^  • The  basic  problem  is  that  while  real 
numbers  are  completely  ordered,  intervals  are  not.  The  result  is  that  no 
entirely  satisfactory  way  of  translating  a comparison  of  real  numbers 
into  a comparison  of  internal  numbers  exists. 

C.  ZEROS  ENTERING  COMPUTED  INTERVALS 

Yet  another  source  of  difficulty  in  interval  arithmetic  stems  from 
the  fact  that  as  a calculation  progresses,  the  widths  of  the  intermediate 
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interval  results  generally  tend  to  become  larger.  In  particular,  to 
guarantee  the  correct  (real)  answer  is  captured  within  the  final  inter- 
val, it  is  necessary  to  round  intermediate  answers  to  machine  representable 
numbers  which  assure  this  inclusion.  While  this  strategy  is  obvious  and 
correct,  the  effect  is  to  expand  interval  width.  In  most  cases  this 
somewhat  conservative  approach  is  defensible  and  does  not  lead  to  prob- 
lems. But  when,  in  the  course  of  a calculation,  the  intervals  become 
null  (contain  the  real  number  aero),  or  overlap,  otherwise  innocent 
operations  may  produce  catastrophic  interval  growth. 

D.  TIMINC  CONSIDERATIONS 

liming  and  space  considerations  can  be  viewed  in  a general  way  as 
well  as  particularizing  them  to  a certain  application  program.  In  this 
paragraph  we  will  comment  in  general  on  these  matters.  In  the  next 
section  of  this  report  more  specific  details  of  timing  will  be  presented 
for  each  benchmark.  The  efficiencies  gained  from  tuning  the  package 
will  be  presented  in  the  section  following  that. 

Because  the  current  version  of  AUGMENT/ INTERVAL  is  designed  for 
generality  and  portability,  there  are  many  inefficiencies  in  the  package. 
Apart  from  the  obvious  increases  in  memory  required  for  interval  over 
real  arithmetic,  there  are  the  time  and  space  costs  of  generating 
function  tables  in  AUGMENT  for  each  use.  Modifying  AUGMENT  to  eliminate 
function  table  generation  with  each  use  was  investigated  and  found  to  be 
a substantial  systems  programming  effort.  However,  in  an  environment 
where  there  is  frequent  production  use  of  AUGMENT,  the  effort  to  do  the 
modifications  would  likely  be  justified.  This  change  has  not  been  made 
in  the  current  system. 
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Another  inefficiency  caused  by  the  desire  for  generality  is  the 
multiple  levels  of  subroutine  calls  involved  in  even  the  simplest 
source  language  statements.  For  example,  the  statement  X = A * B,  for 
A,  B,  X interval  variables,  generates  a total  of  134  subroutine  calls, 
nested  3 deep. 

Because  in  the  H635  version  BPA  is  identified  with  REAL  and 
EXTENDED  with  DOUBLE  PRECISION,  we  were  able  to  reduce  this  sort  of 
inefficiency  by  inserting  replacement  statements  for  subroutine  calls. 
This  work  is  described  in  Section  IV  of  this  report.  Further  effort 
along  the  lines  of  removing  the  nested  subroutine  calls  could  prove 
very  beneficial  to  overall  efficiency  of  the  package. 

Finally,  some  discussion  with  a faculty  colleague  has  concerned 
the  hardware  implementation  of  interval  arithmetic.  There  is  no  doubt 
that  significant  improvements  in  run  times  can  be  achieved  by  replacing 
the  arithmetic  primitives  with  hardware.  However,  dependency  problems 
will  not  be  rectified  by  hardware  implementation. 


VIII.  THE  ENGINEERING  ALGORITHMS 


A.  SPLINE 

A report  made  March  15,  1977  contains  the  results  of  the  benchmark 
run  of  the  subroutine  SPLINE.  These  results  are  also  included  in  this 
final  report  for  the  sake  of  completeness.  As  reported  earlier,  SPLINE 
ran  under  INTERVAL  with  relative  ease.  The  interval  widths  of  the 
answers  varied  from  0.0003  to  0.0012  and  the  intervals  were  very  nearly 
centered  on  the  real  results.  Thus  the  SPLINE  algorithm  seems  to  be 
fairly  stable  when  performed  in  interval  arithmetic  and  as  remarked  in 
Section  II  of  this  report,  this  indicates  stability  of  the  algorithm 
under  real  arithmetic.  The  timing  factor  was  13.4,  that  is,  it  took 
13.4  times  as  long  to  run  under  AUGMENT /INTERVAL  as  it  did  to  run  in 
real  arithmetic.  After  tuning  the  system,  as  described  in  Section  IV 
of  this  report,  the  same  run  was  repeated  with  the  result  that  the 
intervals  were  generally  narrower  and  the  timing  factor  was  reduced  to 
12.2,  a 10%  improvement.  Listings  of  these  runs  are  included  in 
Appendix  B. 

3.  FFT 

Initial  efforts  to  run  the  FFT  benchmark  failed  as  a result  of  an 
error  in  the  program  modifications  written  at  WES.  When  this  problem 
was  finally  ferreted  out  in  June,  the  FFT  results  were  still  not  very 
good.  Closer  examination  revealed  an  error  in  the  INTERVAL  cosine  routine 
sent  to  us  by  MRC.  After  rewriting  INTC0S  we  discovered  that  the  same 
fix  had  been  distributed  by  MRC  last  summer  and  had  been  added  to  our 
working  system.  It  is  not  clear  how  the  correction  was  lost,  nor  whether 
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the  uncorrected  version  was  included  in  the  package  delivered  to  WES  in 


November.  However,  the  recently  delivered  system  contains  the  corrected 
INTCOS. 

Because  of  the  difficulties  mentioned  above,  a smaller  problem  had 
been  constructed  at  K.U.  to  test  FFT.  It  consists  of  an  8-point  approxi- 
mation to  sin  — | over  the  interval  [0,2],  Interval  results  for  this 
test  appeared  at  first  to  have  failed  to  capture  the  real  answer,  that 
is,  the  'correct'  answers  were  not  in  the  intervals  computed.  After 
careful  analysis  of  the  algorithm  and  the  programs,  it  was  determined 
that  the  problem  was  in  the  real  algorithm.  The  intervals  had,  in  fact, 
captured  the  real  arithmetic  results.  But  the  intervals  were  narrower 
than  the  error  in  the  real  number  answers,  and  hence  did  not  extend  far 
enough  to  capture  the  'true'  answers  (see  results  in  Appendix  B) . In 
this  way  the  interval  analysis  did  in  fact  highlight  the  inaccuracies 
of  the  algorithm.  However,  these  results  also  indicate  that  FFT  is  a 
numerically  stable  algorithm  for  this  data.  Going  to  double  precision 
would  probably  improve  accuracy  of  the  real  answers  enough  to  have  the 
intervals  capture  the  'correct'  solutions.  Unfortunately  this  conjecture 
cannot  be  tested  since  double  precision  is  not  available.  But  the  con- 
clusion regarding  the  stability  of  the  FFT  algorithm  remains  valid. 
Listings  and  output  from  the  8-point  runs  is  included  in  Appendix  B. 

There  was  one  significant  change  required  in  the  benchmark  driver. 
The  program  sent  from  WES  included  a calculation  of  the  modulus  of  the 
answer.  This  calculation  being  the  sum  of  squares  of  the  real  part 
(interval)  and  complex-part  (interval)  introduced  a dependency  error  in 
the  form  of  a negative  argument  for  the  square  root  function,  INTSQT. 

To  avoid  this  problem,  the  modulus  calculation  was  taken  out  of  the 
driver. 


47 


c. 


CAUSSE 


The  first  attempts  to  run  GAUSSE  with  the  data  provided  by  WES  all 
resulted  in  aborts  due  to  timer  runouts.  Even  at  a factor  of  100,  that 
is  allowing  the  interval  arithmetic  run  100  times  as  much  time  as  the 
real  arithmetic  run  required,  the  program  was  completing  only  about  two 
thirds  of  the  problem.  Since  this  factor  alone  would  deter  using  GAUSSE 
in  an  interval  arithmetic  setting,  it  was  decided  to  investigate  alterna- 
tive means  for  analyzing  the  algorithm.  At  this  point  we  were  unable  to 
distinguish  the  causes  of  our  problems,  although  it  appeared  that  all  of 
the  difficulties  described  in  Section  II  were  contributing.  The  data 
sent  from  WES  was  a randomly  generated  set  of  coefficients  and  right-hand 
constants  for  a 100  x 100  linear  system,  which  we  discovered  later  was 
incomplete  (the  last  card  had  been  lost  from  the  deck).  When  we  reduced 
the  problem  to  a 20  x 20  array  of  randomly  generated  input,  the  program 
ran  to  completion  but  produced  useless  results  as  reported  in  the  Phase 
111  report  in  March. 

To  gain  better  insight  into  the  nature  of  the  problems  using  INTERVAL 
with  this  algorithm,  we  decided  to  generate  coefficients  for  smaller 
arrays  which  had  properties  known  to  produce  ’good’  results  for  real 
arithmetic  GAUSS  elimination  with  partial  pivoting,  and  then  to  vary 
the  data  toward  less  well  conditioned  problems.  More  than  fifty  runs 
were  made  using  coefficient  arrays  generated  as  follows: 

Upper  left-hand  element  is  called  I0NE0NE 
Diagonal  elements  Increase  by  1 going  down  the  diagonal 
First  super  diagonal  contains  -10 
Second  super  diagonal  contains  -5 
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All  other  elements  are  -1 

Solution  vector  is  all  +1 

By  varying  IONEONE  between  0 and  99900,  we  were  able  to  study  the  behavior 
of  the  algorithm  while  controlling  the  condition  of  the  problem.  Several 
patterns  emerged  quite  clearly. 

(1)  For  well  conditioned  matrices,  the  interval  widths  are  small 

and  actually  improve  (become  narrower)  from  the  last  to  the  first  computed 
result.  There  is  no  pivoting.  See  run  with  IONEONE  = 250  in  Appendix  B, 
for  example. 

(2)  For  less  well  conditioned  matrices,  pivoting  usually  takes 
place  changing  the  order  of  solution,  so  it  is  not  easy  to  follow  the 
accuracy  of  successive  solutions.  The  pattern  that  consistently  emerges 
is  that  the  intervals  become  larger  starting  with  tiie  last  variable  and 
progressing  toward  the  first.  However,  the  intervals  decrease  in  width, 
starting  five  or  six  values  before  the  variable  for  which  pivoting  takes 
place;  and  then,  following  the  pivoted  variable,  increase  exponentially 
(see  example  with  IONEONE  = 2,  size  = 40). 

(3)  For  randomly  generated  data,  the  interval  results  are  worse. 
ror  exa:  le , one  20  x 20  case  produced  intervals  of  virtually  (-<*>,«>) 
for  all  variables,  while  the  best  random  20  x 20  case  managed  only  3 or 
4 place  accuracy. 

Regarding  timing,  well  conditioned  systems  ran  under  INTERVAL  in 
about  16  - 18  times  the  real  run  time  using  the  untuned  system  (reported 
in  March).  With  the  tuned  version  of  the  svstem,  this  factor  was  reduced 
to  about  13.  There  was  significant  increase  in  time  for  INTERVAL  runs 
using  randomly  generated  data.  For  example,  the  20  x 20  random  case  took  15 
times  the  real  case  run  time. 
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IX.  TUNING  THE  SYSTEM 


A.  THE  AUGMENT  PRIMITIVES 

After  studying  the  documentation,  it  was  decided  to  rewrite  the 
eight  machine  dependent  primitives  PACK,  CCODE,  MOVHOL,  NUMIN,  ORDER, 
STRCHR,  STRWDS  and  STRLNC.  These  primitives  had  been  written  for  the 
G635  at  the  MRC  and  were  not  in  the  most  efficient  form.  In  particular, 
use  of  such  functions  as  FLD  was  frequent,  and  produced  grossly 
inefficient  object  code  on  the  G635.  All  eight  primitives  were  rewritten 
in  assembly  language  and  i-  :allec  in  the  package.  See  Appendix  C. 

3.  LINKING 

In  dealing  with  the  space-ef t i c ienc \ balancing  problem,  the  size 
of  the  program  AUGMENT  dictates  the  need  for  overlaying  it  to  reduce 
costs  and  space  requirements.  This  i . sential  In  : small  memory 
installation  but  is  desirable  in  any  a-  version  of  the  system 

delivered  with  this  final  report  is  appro  . at>  • linked  for  overlaying. 

C.  PROVISION  FOR  MISSING  FORTRAN  BKaFi  SUBROUTINES 

The  AUGMENT/ INTERVAL  package  requires  the  availability  of  a specific 
set  of  subroutines  in  the  host  FORTRAN  library.  Honeywell's  FORTRAN 
' ibr;-rv  'or  the  bOO  ' ine  dees  not  include  the  following  required  routines 
in  a readily  usable  form: 

SPA  (REAL)  EXTENDED  (DOUBLE  PRECISION) 

TAN  DTAN  DL0G2*  QTANH 

CBRT  DARSIN  DSINH  DCERT 

CARCOS  DCOSH  DEXP2** 

* Rename  of  DLOG 
♦♦Rename  of  DEXP 


50 


In  order  to  expedite  implementation,  dummy  routines  were  inserted 
for  these  eleven  functions.  Subsequently,  Honeywell  released  a new 
FORTRAN  library  on  their  6600  line  which  included  all  of  the  above 
routines.  However,  the  possibility  exists  that  these  newly  released 
routines  employ  instructions  from  E.I.S.  which  are  not  available  on 
the  WES  computer. 

D.  MACHINE  CONSTANT  FRACBD 

FRACBD  is  a machine  dependent  constant  used  in  various  INTERVAL 
routines.  It  is  intended  to  provide  a direct  means  for  truncating  the 
decimal  part  of  a real  number.  The  value  assigned  to  FRACBD  in  the 
first  implementation  was  close  to  but  not  exactly  equal  the  number 
expected.  As  a result,  there  were  occassional  unexplained  inaccuracies 
in  the  computations.  The  correct  value  has  been  inserted  in  the  current 
version  of  the  package.  (See  Appendix  A). 

E.  MODIFICATION  OF  INTMUL 

In  the  subroutine  INTMUL,  as  a result  of  some  questionable  coding, 
the  variable  TEMP  could  have  been  used  in  a logical  comparison  before 
it  was  assigned  a value.  This  occurred  in  the  statement 

IF (CASE  .NE.  5 .OR.  BPATMP(l)  .LE.  TEMP)  GOTO  110 

after  the  line  labelled  104  and  again  in  the  statement 

IF (CASE  .NE.  5 .OR.  B?ATMP(1)  .CE.  TEMP)  GOTO  120 

after  the  statement  labelled  112.  In  both  cases,  if  TEMP  had  not  been 
assigned,  then  CASE  was  not  equal  to  5,  and  so  the  right  side  of  the 


1 


condition  did  not  have  to  be  checked.  However,  the  Honeywell  FORTRAN 
compiler  does  execute  the  code  which  compares  BPATMP(l)  and  TEMP.  Hence, 
the  contents  of  TEMP  could  cause  errors.  By  assigning  zero  to  TEMP 
upon  entry  to  1NTMUL,  the  possibility  of  error  is  eliminated. 

F.  REMOVING  UNNECESSARY  CALLS  TO  BPA  SUBROUTINES 

Two  changes  were  made,  to  the  description  deck  to  increase  the 
efficiency  of  the  INTERVAL  package. 

First,  all  calls  to  BPaSTR  were  eliminated  by  removing  the  line 

SERVICE  COPY  (STR) 

from  the  BPA  section  of  the  deck.  The  current  version  of  INTERVAL 
equates  the  typos  REAL  and  BPA,  so  A = B is  used  now  instead  of  CALL 
BPASTR(B.A) . 

Second,  all  calls  to  BPALT,  BPALE,  BPAEQ,  BPANE , BPAGE,  and  BPAGT 
have  been  eliminated.  Once  again,  since  REAL  and  BPA  are  the  same, 
the  output  source  from  AUGMENT  uses  the  corresponding  FORTRAN  logical 
operators. 

This  was  accomplished  by  adding  an  ENVIRONMENT  section  to  the 
description  deck,  and  the  FORTRAN  logical  operators  were  defined  for 
BPA  operands. 

To  illustrate  the  effect  of  these  changes,  for  the  statement 
X = A * B,  noted  in  Section  II  D,  the  number  of  subroutine  calls  is 
reduced  from  134  to  46  and  “the  nesting  is  only  2 deep. 
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CODING  OF  THE  PRIMITIVES 
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This  section  contains  a description  of  the  algorithms  developed 
for  BPA  arithmetic  and  for  brounding,  together  with  explanations  of 
the  design  decisions  made. 

Initially  these  algorithms  were  written  to  incorporate  special 
handling  for  the  asymmetric  numbers  described  below.  However,  after 
having  obtained  a working  version  which  allowed  the  special  cases  to 
be  treated  correctly  as  far  as  INTERVAL  was  concerned,  it  was  decided 
these  special  cases  caused  inefficiencies  and  difficulties  of  under- 
standing beyond  their  usefulness.  For  this  reason  a change  was  made 
to  the  algorithms  which  eliminated  the  special  cases  as  either  inputs 
to  or  results  from  the  arithmetic  routines  (there  is  one  exception, 
which  will  be  explained  later) . To  insure  accuracy  during  the 

brounding  of  the  special  cases,  the  original  routine  required  the  input 
to  have  the  correct  sign.  This  requirement  is  no  longer  necessary, 
and  so  the  current  brounding  routine  requires  non-negative  input,  and 
the  sign  correction  is  made  at  the  end. 

Normalized,  floating-point,  two's  complement  numbers  on  the 
Honeywell  635  are  not  symmetric  about  zero.  For  single  precision, 
normalized  numbers  the  following  ranges  are  given: 

negative  -(1  + 2~26) (2_129)  ? N > -2129 

positive  2"129_5  N < (1  - 2_27)(2127). 

“128 

Furthermore,  zero  is  represented  as  (0) ( 2 l). 

Two  main  problems  are  caused  by  this  asymmetry:  the  absolute 
value  of  the  negative  number  with  the  largest  magnitude  cannot  be 
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represented,  and  the  smallest,  normalized,  positive  number  has  no 
normalized  complement  (although  it  does  have  an  unnormalized  complement). 
As  mentioned  above,  these  two  numbers  are  not  allowed  as  ii.puts  to  the 
arithmetic  routines.  Thus  every  valid  argument  has  an  exact  complement 
in  the  INTERVAL  package. 

No  arithmetic  operation  takes  place  if  either  ! these  two  nuabers 
is  found.  When  the  most  negative  number  is  in  argument  l 'he  it  i » « • ■ . 

the  infinity  or  the  overflow  Indicat  >r  is  set,  am:  t;u  « ga  r..  ,n  . r 

with  the  largest  magnitude  allowed  1 reiur  . : is  t • msw>-i 
the  smallest  positive  number  is  used  as  an  argument.  ' 

Indicator  Is  set,  and  either  zero  or  the  sm.i  - • 

is  allowed  is  returned  as  the  answer,  depending  -n  hi  i tin. 

The  single  exception  to  this  treatment  oi  ■ ur - in  : i ■.  >■ 

in  which  the  unchanged  dividend  is  ret  .rue-  i-  ... 

zero  divisor  indicator  is  set  . 

Understanding  the  following  routines  requires  ■ in  w.t  .. 
how  the  floating-point  numbers  are  represented. 

A positive,  single  precision,  normalized,  floating-point  number 
has  an  8-bit  exponent  and  a 28-bit  mantissa.  The  exponent  port  holds 
a two's  complement  Integer,  E,  and  the  mantissa  part  holds  a fractional, 
two's  complement  number,  M.  The  relation  for  a floating-point  number 
Z is  this: 

Z = M * 2E. 


Furthermore,  the  mantissa  must  have  a zero  in  the  first  bit  (the  sign 

bit),  and  a 1 in  the  second  bit.  Thus  if  M is  positive,  then 
-1  -27 

2 M _<  1 - 2 . If  we  catenate  non-zero  bits  to  the  right  end 
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of  a positive  mantissa,  then  the  magnitude  of  the  number  represented  is 
increased;  if  we  truncate  non-zero  bits  from  the  right  end  of  a 
positive  mantissa,  then  the  magnitude  is  decreased. 


Double  precision  numbers  also  have  8-bit  exponents,  but  their 
mantissas  are  expanded  to  64  bits.  The  arithmetic  and  brounding 
routines  use  the  full  AQ  register  to  manipulate  the  M(A) , and  so  a 
single  precision  mantissa  occupies  M(A)q  _ 27  with  MCA^g  72  being 
initially  set  to  zero.  The  U register  is  simulated  by  a word  pair  in 
memory. 

Given  the  elimination  of  the  normalized  numbers  without 
(normalized)  complements,  the  range  of  numbers  allowable  in  the 
INTERVAL  package  is  the  following: 

negative  -(1  + 2~26) (2~129)  >N>-(1-  2~27)(2]27) 

|j  positive  (1  + 2-26) (2~129)  £ N < (1  - 2~27) (2127) . 

Four  machine  dependent  constants  are  used  in  the  following  routines. 


MINNEG  is  the  negative  number  with  the  largest  magnitude. 

MINPOS  is  the  smallest  positive  number. 

POSBND  is  the  smallest  positive  number  allowed  in  the 

INTERVAL  package. 

MAXPOS  is  the  largest,  single  precision,  positive  number. 

In  addition,  two  addition  constants  are  used  by  the  brounding  routine, 
and  their  mantissas  are  given  below. 


t 


1 

1 


ROUND 


AUGMNT 


00  . . 

. 00 

100  . 

...  00 

0 

27 

28 

63 

00  . . 

. 01 

00  . 

...  00 

0 

27 

28 

63 
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When  one  of  these  addition  constants  is  used,  its  exponent  is  set  to 
the  exponent  of  the  number  to  which  it  is  being  added.  The  mantissas 
of  ROUND  and  AUGMNT  are  not  normalized. 

The  arithmetic  routines  BPACEB,  BPASUB,  BPAADD , BPAMUL,  and 
BPADIV , and  the  brounding  routine  BROUND  are  given  below.  Additional 
comments  are  given  for  the  numbered  lines. 


BPACEB  : <<Convert  the  double  precision  A1  to  single 
<<precision  and  bround  it. 

<<RES  +-  A1  >> 

clear  all  indicators  ; 

A •*-  A1  ; 
case 

: A = MINNEG  : <<Overflow.>> 

SC  1 ; goto  .EO.  ; 

: A = MINPOS  : <<Underf low. >> 
goto  ,EU . ; 
endcase  ; 
if  A < 0 

then  A ■<-  -A  ; SC  1 ; 
fi  ; 

goto  BROUND  ; <<End  of  BPACEB>> 


BPASUB  : <<Compute  RES  A1  - A2.>> 

clear  all  indicators  ; 

U A2  ; 
if  U = 0 

then  <<The  answer  is  Al.>> 

A * A1  ; 
case 

: A.=  MINNEG  : <<Overf low.>> 
SC  *■  1 ; goto  .EO.  ; 

: A = MINPOS  : <<Underf low.>> 
goto  .EU.  ; 

endcase  ; 

RES  ■*-  A ; return  ; 


case 

: U = MINNEG  : «Overflow.» 

SC  1 ; goto  .EO.  ; 

: U = MINPOS  : <<Underf low. >> 
A -e  U ; goto  .EU.  ; 

endcase  ; 
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U «-  -U  ; 

goto  SUBIN  ; <<BPAADD  completes  the  subtraction. >> 


BPAADD  : «Compute  RES  ■*-  A1  + A2.» 


clear  all  indicators  ; 

U + A2  ; 
if  U = 0 

then  <<A1  is  the  answer. >> 


A = MINNEG  : «Overflow.» 

SC  1 ; goto  .EO.  ; 

A = MINPOS  : «Underf low.>> 
goto  .EU.  ; 


return 


MINNEG 
SC  <-  1 
MINPOS 
A+ll 


<<Overf low. >> 
goto  ,E0.  ; 

<<Underflow.>> 
goto  .EU.  ; 


endcase 


SUBIN  : <<This  is  the  entry  point  for  the  subtraction  routine.  >> 
A •*-  A1  ; 
if  A = 0 

then  <<U  holds  the  exact  answer. >> 

RES  *■  U ; return  ; 


A = MINNEG  : <<Overf low. >> 
SC  *■  1 ; goto  .EO.  ; 

A = MINPOS  : <<Underf low.>> 
goto  .EU.  ; 


if  [U|  < | A | 

then  U ■*-*•  A T 

fi  ; 

ijf  A < 0 

then  <<Make  A positive .'>^ 

A <-  -A  ; U •*-  -U  ; SC  v 1 ; 

f i ; 

TEMP  +•  E(U)  - E(A)  ; <<Compute  shift  taunt. >> 


if  TEMP  > 0 «No  shift  if  TEMP  _<  0.>> 
then 

E(A)  - E(U)  ; 

shift  M(A)  right  TEMP  bits  ; <<Line  up  the  mantissas. >> 

fi  ; 

A <-  A + U ; 
if  A = 0 

then  <<The  answer  must  be  exact. >> 

RES  A ; return  ; 

li  ; 

if  RI  *=  1 

then  <<We  lost  digits  in  the  mantissa,  and  the 
<<number  in  A is  too  small,  regardless  of 
<<whether  it  is  positive  or  negative.  We 
<<set  bit  35  to  1;  this  bit  is  in  the  extra 
<<precision  word,  and  will  not  affect  the 
<<rounding  option.  >> 

M(A)35  <-  1 ; 

fi  ; 

if  A < 0 

then  A *■  -A  ; SC  -e  -SC  ; 

fi  ; 

goto  BROUND  ; «End  of  BPAADD» 


: 


Additional  comments  for  BPAADD: 

Lines 

1 At  the  present  time  (1976-1977)  the  floating  point  compare 
magnitude  instruction  DFCMG  does  not  work  on  the  Honeywell 
635;  consequently,  we  program  around  it  by  making  both 
arguments  positive,  and  then  comparing  them.  The  problem 
with  DFCMG  was  not  known  by  Honeywell  at  the  time  it  was 
reported  to  them. 

2-4  The  exponent  of  the  absolute  value  of  an  exact  power  of  two 
is  alwcys  one  greater  than  the  exponent  of  its  complement; 
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consequently,  E(U)  - E(A)  will  be  negative  if  we  are  adding 
a power  of  two  and  its  complement.  If  the  difference  in  the 
exponents  is  negative  or  zero,  then  no  shifting  need  be  done. 
Otherwise,  we  need  to  live  up  the  mantissas  for  the  addition  and 
change  exponent  of  A accordingly.  The  residue  indicator  may 
be  set  during  the  shift. 

5 If  the  computed  result  is  MINPOS,  then  this  negation  will  cause 

the  exponent  underflow  indicator  of  the  machine  to  be  set, 
and  then  it  is  treated  as  an  underflow-by-one  in  BROUND.  The 
result  can  not  be  MINNEG. 

BPAMl'L  : <<Compute  RES  A1  * A2.>> 
clear  all  indicators  ; 

U +•  A2  ; 

if  U = 0 oir  A1  = 0 

then  RES  •*-  0 ; return  ; 

fi  ; 
case 

: U = MINNEG  : <<Overflow.>> 

SC  *■  1 ; goto  .EO.  ; 

: U = MINPOS  : «Underflow.» 

A U ; goto  .EU.  ; 

endcase  ; 
if  U < 0 

then  <<Make  it  positive. >> 

U + -U  ; SC  + 1 ; 

fi  ; 

A A1  ; 
case 
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: A = MINNEG  : <<Overflow.>> 

SC  *■  1 ; goto  .EO.  ; 

: A =•  MINPOS  : «Under£low.» 

SC  +■  0 ; goto  .EU.  ; 

endcase  ; 

If  A < 0 

then  <<Make  it  positive. >> 

A «-  -A  ; SC  -<-  -SC  ; 

fi  J 

1 A +-  A * U ; 

goto  BROUND  ; <<End  of  BPAMUL» 

Additional  comments  for  BPAMUL: 

Lines 

1 Single  precision  mantissas  have  28  bits  (27  bits  plus  a sign), 
and  so  the  product  can  be  held  in  55  bits  ( a sign  bit  plus 
54  for  the  mantissa).  FMP  computes  the  full  AQ  (72  bits) 
accurately,  and  so  we  use  the  machine's  floating  multiply 
instruction. 

BPADIV  : <<Compute  RES  ■*-  Al/A2.>> 
clear  all  indicators  ; 

0 + A2  ; 

if  U = 0 

I then  <<Division  by  zero.>> 

BPAFLT  + 4 ; «DZ  *■  1» 

IRES  A1  ; <<No  checking  here.>> 

return  ; 

I | 

I 
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fi  ; 

if  A1  - Q 

-rhen  RES  *-  0 ; return  ; 

fi  ; 
case 

: V = MINNEG  : «Overflow.» 
SC  ■*-  1 ; goto  .EO.  ; 

: U = MINPOS  : <<Underf low. >> 
A •<-  U ; goto  .EU.  ; 

endcase  ; 
if  U < 0 

f <<Make  it  positive. >> 

U +■  -U  ; SC  «-  1 ; 

fi  ; 

A «-  A1  ; 
case 

: A = MINNEG  : <<Overflow.>> 
SC  *■  1 ; goto  .EO.  ; 

: A = MINPOS  ; «Underflow.» 
sc  +-  0 ; goto  .EU.  ; 

endcase  ; 
if  A < Q 

— <<:Make  it  positive. >> 

A <-  -A  ; SC  «-  -SC  ; 

fi  ; 

V 0 ; 

f°r  j 35  down  to  1 do 
|M(A)|  _>  | M(U)  J 

then 


. — 


2 

3 


4 


I 


5 

6 

7 

8 

9 


M(A)  +-  M(A)  - M(U)  ; 

V V + 1 ; 

fi  ; 

shift  M(A)  and  V left  by  1 ; «Multiply  by  2.» 
endfor  ; 

shift  M(A)  right  36  bits  ; <<Signal  remainder. >> 

M(A)  V ; <<Move  the  quotient  to  M(A).>> 
if  M(A)  = 1 

then  <<We  need  to  normalize  the  quotient. >> 
shift  M(A)  right  1 bit  ; 

E(A)  *■  E(A)  + 1 ; <<Correct  the  exponent. >> 

fi  ; 

E(A)  +-  E (A)  - E (U)  ; 
if  overflow 

then  <<We  may  have  exponent  overflow  or  underflow.  >> 
if  E(A)q  = 1 

then  <<Overflow. >> 
goto  ,E0.  ; 
else  <<Underflow. >> 
goto  .EU.  ; 

fi  ; 

fi  ; 

if  RI  = 1 

then  <<Set  a'low  order  bit  in  M(A)  to  1.  » 

M(A)35  -e  1 ; 

11  ; 

goto  BROUND  ; «End  of  BPADIV.>> 


| 
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Additional  comments  for  BPADIV: 


Lines 

— 


1-4 


Yohe's  algorithm  for  division  in  the  base  6 number  system 
includes  an  extra  loop  at  this  point  to  subtract  M(U) 
from  M(A)  until  M(A)  < M(U) . The  normalised  form  of 
floating-point  numbers  in  the  Honeywell  635  makes  this 
extra  loop  unnecessary;  the  M(A)  must  always  be  less  than 
2 * M(U),  and  50  the  subtraction  in  line  2 insures  that 
M(A)  < M(U) . 


r 

Li 


5-6  If  the  subtraction  in  line  5 produces  a number  either  too 

large  or  too  small  to  be  represented  in  the  8-bit  exponent 
of  the  635,  then  the  overflow  indicator  of  the  machine  is 
set. 

7-9  To  determine  whether  the  fault  was  caused  by  exponent 

overflow  or  underflow,  we  inspect  the  exponent  computed. 

If  E(A)  is  too  big,  then  the  machine  uses  the  sign  bit  to 
gain  the  magnitude  it  needs.  Thus  exponent  overflow  sets 
E(A)q  to  1.  If  E(A)  is  too  small,  then  E(A)q  becomes  zero 
when  the  subtraction  is  performed. 


BROUND  : <<Bround  the  number  in  A,  and  return  the  answer  in  RES. 
<<The  brounding  option  is  held  in  OPTION.  The  fault 
<<value  is  returned  in  BPAFLT,  and  the  following 
<<values  are  used: 


<< 

0 

operation  successful 

<< 

1 

exponent  overflow 

<< 

2 

infinity 

<< 

3 

exponent  underflow 

<< 

4 

zero  divisor 

— 
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if  OPTION  i {'U',  ’L',  'T',  'R',  'A'} 

then  <<Use  R if  OPTION  is  invalid. >> 


OPTION  ’R'  ; 

11  ; 

if  exponent  overflow  or  A > MAXPOS 
then  goto  . EO . ; 

ft  ; 

if  exponent  underflow  or  A < POSBND 
then  goto  . EU . ; 

11  ; 

<<The  correction  for  the  residue  indicator  has  been  made 
<<already,  and  is  held  in  M(A>3g.  >> 

case 

1 : [(OPTION  = 'U'  and  SC  = 0)  or 

2 (OPTION  = 'L'  and  SC  = 1)  or 

3 OPTION  = 'A']  and  M(A)9g_63  + 0 : 

RES  - A + AUGMNT  ; 
if  exponent  overflow 

then  <<Set  the  infinity  flag  and  load  MAXPOS. » 
3PAFLT  •*-  2 ; 

RES  MAXPOS  ; 


11  ; 


(OPTION 


(OPTION 


OPTION 


then  RES  *-  -RES 


return 


OPTION 


RES  ■*-  A + ROUND 


if  exponent  overflow 


then  <<set  the  infinity  flag  and  load  MAXPOS>> 


EPAFLT  *■  2 


RES  *•  MAXPOS 


<<Overflow  faults  are  handled  here.  The  input 
<<number  Is  irrelevant. 


(OPTION 


(OPTION 


OPTION 


<<The  brounding  option  implied  truncation, 

<<and  so  use  the  fault  value  for  exponent 
<<overflow.  >> 

BPAFLT  «-  1 ; 

else  <<Augmentation  was  implied,  and  so  we  had 

«an  infinity  fault.  » 

BPAFLT  +-  2 ; 


endcase  ; 

RF.S  *■  MAXPOS  ; 
if  SC  - 1 

then  RES  *■  -RES  ; 

11  ; 

return  ; 

.EU.  : <<\jnderflow  faults  are  handled  here.  Special 
<<action  is  taken  if  the  R option  is  specified 
<<and  we  have  exponent  underflow  by  one.  >> 

BPAFLT  •*-  3 ; <<exponent  underflow>> 

case 


: (OPTION  = 'U'  and  SC  = 0)  or 
(OPTION  = ’L’  and  SC  = 1)  or 
OPTION  = 'A'  : 

<<  The  brounding  option  implied 
<<  augmentation,  and  so  we  use  the 
<<  non-zero  number  with  the  smallest 
<<  magnitude  allowed.  >> 

RES  <-  POSBND  ; 

: OPTION  t ’R’  : 

<<Truncation  was  implied;  so  we 

<<use  zero.  >> 


RES  +•  0 ; 
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else  <<OPTION  = 'R'» 


7 J^f  A POSBND  o£  exponent  underflow  by  1 

then  RES  •*-  POSBND  ; 
else  RES  ■*-  0 ; 

fi  ; 

endcase  ; 

if  SC  = 1 

then  RES  •* — RES  ; 

fl  ; 

return  ; «End  of  BROUND>> 

Additional  Comments  for  BROUND: 

Lines 

1-3  Augmentation  is  implied  in  3 cases:  (1)  OPTION  = 'U' 

and  SC  = 0 means  the  true  answer  is  positive,  and  we 
want  the  least  upper  bound;  (2)  OPTION  = 'L'  and  SC  = 1 
means  the  true  answer  is  negative  and  we  want  the 
greatest  lower  bound,  thus  we  must  augment  the  positive 
number  to  increase  the  magnitude;  (3)  OPTION  = 'A'. 

If  !1(a)28-63  = ® then  the  result  is  exact. 

4-6  Truncation  is  implied  in  4 cases:  (1)  OPTION  = ’U' 

and  SC  = 1 ijeans  the  true  result  is  negative,  and  we 
want  the  least  upper  bound,  thus  we  truncate  the  positive 
result  to  decrease  the  magnitude;  (2)  OPTION  = 'L'  and 
SC  = 0 means  the  true  result  is  positive,  and  we  want 
the  greatest  lower  bound;  (3)  OPTION  - 'T';  (4)  M(A)^Q  ,,  = 0 
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Lines 


means  the  result  can  be  expressed  exactly  in  single 
precision,  and  no  brounding  option  will  affect  it, 
thus  the  extra  precision  bits  are  effectively  truncated. 


4-6 


i 


7 Exponent  underf low-by-one  can  occur  two  ways:  by 

computation  and  by  input  of  MINPOS.  All  numbers 
smaller  than  POSBND  and  greater  than  or  equal  to  MINPOS 
are  considered  to  be  cases  of  underflow  by  one.  If  we 
actually  cause  a machine  exponent  underflow  during 
computation,  the  computed  exponent  is  127,  or  the 
largest  exponent  possible. 
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X.  DESCRIPTION  OF  THE  TAPE  SENT  TO  WES 

The  tape  sent  to  WES  is  a multifile  tape,  written  at  800  bpi,  7 track. 
The  fifteen  files  contained  on  the  tape  are: 

File  1 Arithmetic  and  brounding  routines  for  INTERVAL 
File  2 Machine  dependent  primitives  for  AUGMENT 
File  3 SESOL  and  BANSOL 

File  4 INTERVAL  source  library,  including  corrected  1NTCOS  and  the 
modifications  to  the  AUGMENT  description  deck  to  remove  the 
superfluous  subroutine  calls 
File  5 FFT 

File  6 SPLINE (INTERVAL  form) 

File  7 CAUSSINT  (old  version  of  CAUSE) 

File  8 GAUSS  (INTERVAL  form) 

File  9 SESOL 

File  10  FFT  (original) 

File  11  FFT  (no  complex) 

File  12  FFT  (INTERVAL) 

File  13  FFT  (8  point  real) 

File  14  FFT  (8  point  INTERVAL) 

File  15  GAUSS  (a  third  version) 
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XI.  CONCLUSIONS  AND  RECOMMENDATIONS 


In  ray  opinion,  the  most  significant  single  fact  that  has  emerged 
from  this  project  is  that  interval  arithmetic  has  limited  value  as  a 
tool  for  analyzing  real  algorithms.  The  limitation  is  specifically 
dependency.  As  pointed  out  by  our  results  and,  in  a more  general  way, 
by  Ris  [2],  the  fact  that  interval  arithmetic  gives  intervals  which 
are  too  wide  does  not  tell  you,  necessarily,  that  the  real  algorithm 
is  bad  nor  that  the  real  algorithm  is  necessarily  sensitive  to  the  real 
data  used.  In  particular,  it  is  well  known  that  Gauss  elimination  with 
partial  pivoting  is  a stable  solution  technique  in  a real  number  setting 
provided  one  is  careful  £ bout  conditioning  and  error  control.  What  is 
not  so  well  known  is  that  the  same  Gauss  method  is  subject  to  severe 
dependency  problems  which  can  cause  as  much  as  fourfold  increases  in 
interval  widths  at  each  stage  of  the  procedure  when  employing  interval 
arithmetic.  The  point  is  that  interval  analysis  yields  information  about 
interval  algorithms  in  all  cases;  but  interval  analysis  yields  information 
about  the  real  version  of  an  interval  algorithm  only  when  the  interval 
algorithm  is  'good'.  One  needs  a way  to  break  down  the  sources  of 
increased  width  when  the  computed  intervals  are  too  wide  in  order  to  be 
able  to  say  any tiling  about  the  corresponding  real  algorithm.  Even  then, 
the  conclusion  may  be  that  nothing  can  be  inferred  from  the  interval 
analysis.  Both  Ris  [2]  and  Hansen  [1]  attack  this  problem,  but  from 
different  points  of  view.  Hansen  proposes  a generalized  interval  arith- 
metic (g.i.a.)  to  reduce  the  effect  of  dependency.  By  standardizing 
the  form  of  interval  representation  and  redefining  the  arithmetic 
operators,  he  is  able  to  reduce  the  inherent  lack  of  sharpness  (caused 
by  dependency)  to  a second  order  effect.  Ris  tackles  the  problem  in  a 
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different  way,  by  defining  predicates  and  functions  on  intervals  which 
have  predictable  propagation  properties  in  interval  arithmetic  operations. 
He  then  applies  these  ideas  to  analyze  interval  algorithms,  highlighting 
those  characteristics  which  might  be  used  to  predict  the  success  or  failure 
of  the  algorithm.  Hansen's  approach  is  more  in  keeping  with  the  current 
project  in  that  it  reduces  the  factor  (dependency)  that  limits  the  value 
of  interval  analysis  of  real  algorithms.  Ris'  approach  is  designed  to 
develop  new  algorithms  for  which  interval  arithmetic  is  well  suited. 

Of  course,  once  a 'good'  interval  algorithm  is  produced,  it  will  be  a 

/ 

/ 

'good'  real  algorithm.  / 

/ 

/ 

Getting  more  specific,  the  results  of  the  benchmark  runs  indicate  j 
the  following:  / 

SPLINE  is  reasonably  stable  for  the  data  run  / 

FFT  is  very  stable  but  suffers  from  loss  of  accuracy  in  the/ 

calculations  for  the  4096  point  case.  / 

/ 

GAUSSE  is  (predictably)  stable  for  well  conditioned  matrices 
and  becomes  less  so  as  the  condition  of  the  matrix 
deteriorates.  For  randomly  generated  data,  the 
algorithm  is  less  efficient  and  less  accurate. 

Pivoting  greatly  increases  widths  due  to  dependency. 

Recommendations: 

/ 

In  view  of  the  number  and  range  of  difficulties  encountered  in 
attempting  to  use  this  AUGMENT/ INTERVAL  package  for  analyzing  algorithms, 
it  clearly  should  be  treated  as  an  experimental  tool,  and  not  considered 
for  routine  or  production  use  in  its  present  form. 

The  magnitude  and  complexity  of  many  of  the  problems  being  solved 
at  WES  makes  interval  analysis  of  the  solution  algorithms  and  sensitivity 
analysis  of  the  data  extremely  difficult  to  accomplish  with  very  high 
reliability  using  the  current  package.  A major  effort  to  employ  interval 


arithmetic  might  better  be  directed  toward  development  of  new  algorithms 
based  on  interval  concepts  (as  implied  in  Ris'  work)  than  to  attempt 
analysis  of  the  real  algorithms  currently  in  use.  Alternatively, 
redefining  the  arithmetic  operations  and  interval  representation  (as 
Hansen  suggests)  might  convert  INTERVAL  into  a more  practical  tool  for 
analysis  of  the  real  algorithms. 
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The  following  are  additional  computer  listings  and  runs  available  from 

the  Automatic  Data  Processing  Center,  U.  S.  Army  Engineer  Waterways  Experiment 

Station,  P.  0.  Box  631,  Vicksburg,  Miss.  39180: 

Interval  Primitives  written  at  the  University  of  Kansas 
Augment  Primitives  written  at  the  University  of  Kansas 

Test  programs  and  their  results: 

Summation  of  first  129  terms  of  (1/X)**(I-1) 

Roots  of  a quadratic  equation 
635  rounding  errors  in  addition 
Test  mathematical  functions  using  INTRAP 
Gaussian  elimination  routine- interval  extension 
Table  of  factorials  and  their  natural  logarithms 
SPLINE  (real  and  interval) 

FFT  (8-point,  real  arc  interval) 

FFT  (4096-point,  real  and  interval) 

Gaussian  (20X20,  real  and  interval) 


k. 
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APPENDIX  A:  MACHINE  DEPENDENT 
CONSTANTS  FOR  'INTERVAL  II' 


1 


positive  BPA  number  such  that  the  low-order  digit  immediately 
the  radix  point. 
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■ 
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APPENDIX  B:  PRIMITIVES  WRITTEN  AT  THE 
UNIVERSITY  OF  KANSAS 


% 

'4 


1 


$ 

fihAP 

00000607 

SYMDE  F 

RP  AA  DO/BPASUR/ 

BP  A D 1 V ,BPAMUL,BPAC£b 

C000C60 3 

BLOCK 

BP  AC OM 

UOOOO  o09 

OPTION 

"SS 

1 

THIS  IS  FOR  GROUNDING  OPTION 

00000610 

BPAFC  T 

9 S S 

7 

THIS  JILL  RETURN  ERROR  CODE 

0000061 1 

USE 

0000061 i 

R SSA 

BBSS 

8 

00000613 

FRLK 

0000061 4 

SHIFT 

BSS 

1 

Tr00~ff0  6l  5 

TEST1 

EBSS 

2 

OOOOC  61  6 

TEST? 

BSS 

2 

0000061 7 

0KTA9  CONTROLS  BRANCHING  IN  BROUNO. 

oooaooi s 

UPPER  HALVES  0 F THE  WORDS  HOLD  BRANCH 

CQOOC61 9 

VALUES  FOR  THE  POSITIVE  NUMBERS. 

00000620 

THE  LOWER  HALVES  HOLD  JUMP  VALUES 

0000  C 62 1 

FOR  THE  NEGATiv/E  NUMBERS. 

0 OOCO  0?  2 

OK  TAB 

OCT 

0 

THE  FIRST  ENTRY  IS  E ^ P T Y 

00000623 

OCT 

r 

UC+)  = 00*  U(-)  = 00  IN  DECIMAL 

00000624 

OCT 

OCOO  30000030 

L<+>  = 24,  L(->  = 24  IN  DECIMAL 

00000625 

OCT 

OOCO JOOOOOOC 

T ( ♦ > = 24,  TC-)  = 00  IN  DECIMAL 

00900626 

OCT 

‘ C00D1  3Cr00322~ 

Hi'*)  = t r,  R(-T  = 1g  IN  DECIMAL 

0 003  3 62  7 

OCT 

9C09000C003U 

A ( ♦ ) = 00,  A ( — ) = 24  IN  DECIMAL 

0000062  i 

EOTAB  IS  A TABLE  WHICH  HOLDS  THE  CORRECT 

00000629 

FAULTS  FOR  OVERFLOW  8R0UND 1NG.  VALUES 

00000630 

FOR  POSITIVE  NUM3ERS  ARE  IN  THE  UPPER 

0000U631 

HALVES  OF  THE  WORDS/  AND  ThE  VALUES  FOR 

00000652 

NEGATIVE  NUMBER'S  ARE  TN  THE  TU'wET* 

C 30(3  0 633 

HALVES.  1 IS  OVERFLOW.  2 IS  INFINITY. 

003CQ  634 

TOTAS 

OCT 

0 

first  Entry  is  empty 

COTOO  63  5 

OC  T 

C000020C0001 

U(+)  = INF,  UC-)  = ov 

00000636 

OCT 

OCCOOIOOOOO? 

L (+)  = OV,  LT-)  '=  INT  ' 

CUOT.  63T" 

OCT 

OOOCCTOGOOOI 

T (+ ) = OV,  T (-)  = OV 

0 J00063S 

OCT 

OVOVVTVnW  D2 

R CO  * T ST,  ET-T  = J'fTF 

GOD0  0 6T9- 

OCT 

000002000002 

A(*>  = INF,  A ( - ) = INF 

00000640 

* 

EUTABP  AND  EUYAflN  ARE  TABLES  'WHICH 

0 COO  0 64  1 

# 

HOLD  THE  CORRECT  ANSWERS  FOR  POSITIVE 

00030 j42 

• 

. _ ... 

AND  NEGATIVE  'JNDERFlOV  RESPECTIVELY. 

0C3C0  64  3 

* 

eu  STANDS  for  exponent  UNDERFLOW. 

00000644 

EUTA9P 

OCT 

0 

f irst  TNTirY  rs  ttptv 

ClOOOrj  64  $ 

OCT 

400400000001 

U C ♦ > = POSBNO 

00000646 

OCT 

400000000000 

L C + 0 = ZERO 

00000647 

OC  T 

4C0000000000 

T(»)  = ZERO 

00030648 

OCT 

400000000000 

R (+ ) = ZERO  EXCEPT  FOR  EU  BY  1 

C 0000  64  9 

OCT 

4 004  0000000 1 

A (+)  = POSBNO 

OCOCC  650 

EUTABN 

OCT- 

p 

f iRsT"rNiR y rrTMBfY 

0 0 r 3 j 6 $ i 

OC  T 

400000000000 

U(-;  = ZERO 

00000652 

OCT 

401377777777 

L(-)  * MAX  NE  G 

0000065 3 

OCT 

400000000000 

T(-)  = ZERO 

00000654 

OCT 

400000000000 

RC-)  = ZERO  EXCEPT  FOR  EU  BY  1 

“00000655 

OCT 

401377777777 

AC—)  = MAX  NE  G 

00000656 

ZERO 

EOCT 

4 corotriTtmcroTr 

“ ZERO  U TT*  C?*  *'  -T2R0 — 

u odo  3 6S  y 

OCT 

0 

MANTISSA  OF  DBL  PREC  IS  0.  THIS  IS  ZERO 

00000658 

I 


cmk  dbi 

OC  T 

000000000377 

USED  TO  CHECK  FOR  ONES  IS  DBLE 

00000659 

1 

OC  T 

777777777777 

PRECISION  PART  gF  MANTISSA 

00000  66  J 

M AXPOS 

OCT 

576777777777 

EXP  IS  1?7<L4RGE$T  POSS.)  MANTISSA 

0 OOUO  66 1 

OC  T 

7777 77777777 

Is  LARGEST  POSITIVE  FRACTION. 

00000662 

MI NPOS 

OC  T 

4004  COOOOOCC 

E X°  IS  -12?,  SO  THIS  IS  MIN  POS. 

00000663 

OC  T 

0 

QQQQ0664 

P OS  BN  D 

OC  T 

AC04  ooocoeoi 

SMALLEST  positive  NUMbER  ALLOWED 

00000665 

oc  r 

c 

00000666 

MA  KNEG 

OC  T 

A L 1 3 77777777 

EXP  IS  -123  AND  MANT  IS 

00000667 

oc  r 

777777777777 

SMALLEST  NEG. 

00000663 

M I NNE G 

OCT 

3 7 7P0O0C  OQ00 

EXP  IS  127  AND  MAN  IS  NEG,  NORMALIZED 

00000669 

OCT 

0 

MINNEG  NOT  ALLOWED  AS  ARGUMENT  OF  ANSWER 

00000670 

NEGBND 

OCT 

3 7 7OO0OOOU J 1 

NEGBND  IS  SINGLE  PREC.  NEG  OF  MAXPOS 

00000671 

OCT 

n 

AND  IT  IS  THE  MOST  NEGATIVE  NO.  ALLOWED 

00000672 

UT  A 

OCT 

p 

USED  FOR  IJ,T  AND  A OPTIONS 

C 0000 673 

OC  T 

777777777777 

OOOOC  674 

POSRND 

OCT 

C 

ADO  CONST  FOR  R ♦ 

COOOC675 

OCT 

4000000000UC 

C 0000  6 76 

NEGRND 

OC  T 

0 

ADD  CONST  FOR  R- 

0000 067 7 

OC  T 

377777777777 

00000678 

EUONE 

OC  T 

376000000000 

MAX  EXPONENT,  UNDERFLOW  8Y  1 

0000  0 6 7 6 

RSTA9L 

OC  r 

777,1 777,5777 

,7777,1 7777,37777,77777,1 77777,377777 

00000680 

OCT 

77 77 77,  1 7 777  77, 3 7777  77,  777  7777, 1 77  77 7 7 7, 3 7 777777,  7 777  7 7770 OOC 0681 

OCT 

1 77777777,3  7 7 

777 777, 777 777777,1 7 7 7777 777, 3777777777 

00000632 

OCT 

7777777777,17777777777,37777777777,77777777777 

C 000  0 68  3 

OCT 

1 77777777777, 

377777777777,777777777777 

00000684 

ADDONE 

EOCT 

n 

00000685 

OCT 

oci^oonooooo 

CC000636 

SUBONEEOCT 

0 

00000687 

OCT 

1 

000006:8 

8IGCHK 

OCT 

376777777777 

USED  TO  TEST  FOR  OVERFLOW 

COCOO  66  9 

OC  T 

0 

AGAINST  MAXPOS. 

00000670 

EUOOF  F 

OCT 

74‘777'r 

00600691" 

EOON 

OCT 

C2COOO 

EXPONENT  OVERFLOW  IS  BIT  22 

00000672 

E U ON 

OC  T 

010000 

EXPONENT  UNDERFLOW  IS  BIT  23 

00000693 

OF  OMSK 

OCT 

000000004000 

BIT  24  13  THE  OVERFLOW  MASK  INDICATOR 

C 000  0 694 

TEMP 

ERSS 

? 

00000695 

ARG1 

E 0 S S 

2 

00000696 

ARG2 

BS  S 

2 

00000697 

COUNT 

s 

1 

00000678 

SAVEXP 

BSS 

1 

00000699 

SAFE 

EDSS 

2 

00000700 

R .U. 

PSS 

2 

00000701 

ANS 

BSS 

2 

00000702 

BITS 

oc  T 

1 C. 00 000/? GUO0CC/ 4 03  0 GuG# 10000000# 2000 OOOC, 4 0000000 

CO  00  0 70  3 

OCT 

1 00000000,2 00000000,4 00000000,1000000003, 2000000000 

00000704 

OCT 

AOOOOOCOOC, 1 A 

OOF  OCT'  2 00, 2 '■’00  COCDO  00, 4 C 0000  OCQOO 

COOOC705 

OC  T 

1 OOrCOCOOOOO, 

200000000000,400000000000 

00000706 

CALL 

OSS 

1 

USED  TO  IDENTIFY  CALLING  ROUTINE 

00000  70  7 

TURN 

MACRO 

00000  708 

STQ 

TEMP* 

00000709 

LDQ 

ffl 

00000710 

f 


[ 


ENTER 


ST  I 

TEMP 

OCOUG  71 1 

OR  SQ 

TEMP 

GOO ju  71  2 

L D I 

TEMP  ' 

~cjm  j7i  3 

£ N DM 

TURN 

OOOOC  71  4 

MACRO 

UOOOO  71  5 

SREG 

. RSS  A 

00000  71  6 

ST  I 

• E .L  . . 

C'OCCJ  71  7 

S T X 1 

.E .L  . . 

00000  71  6 

rocnnrm- 

L D A 

OF  OMSK 

WE  "WANT  TO  TURN  OFF  TOE  F A ULT- TRaF5 

ST  I 

T EMP 

GET  A COPY  OF  THE  INDICATORS 

CC00072C 

ORSA 

TEMP 

SET  OVERFLOW  MASK  INDICATOR  OF  COPY 

COOOO  721 

LD  A 

EUOO  F F 

TURN  EXPONENT  OFLO  AND  UNFLO  OFF 

UGOOC  722 

ANSA 

TEMP 

00000723 

LDI 

TEMP 

FAULT  TRAPS  ARE  NOW  DISABLED 

C 000  0 724 

LD  X2 

0/  DU 

CLEAR  ERROR  REGS. 

0 GOO 0 72  5 

LDX4 

0 / DU 

00000  726 

L D X 6 

0,  DU 

CLEAR  SC 

00000727 

LDX7 

p,du 

C Lr  A R R I 

00000728 

F A x r* 

2,1  * 

GET  FIRST  ARG 

00000729 

I FE 

*1/1 ,21 

IF  A DD, SU8/MUL/0R  DIV,  DO  NEXT  21  LINES 

COOOO  730 

FLO 

■'0/0  " 

u 0 3 0 u 7 3 1 

FCMP 

M I NN E G 

“INNEG  IS  NOT  ALLOWED  AS  AN  ARGUMENT 

00000732 

TN2 

* + 2 

OOOOC 73  3 

TR  A 

8R0UND 

AND  GO  TO  SROUND  PART  .EO. 

00000754 

F C MP 

MI NPOS 

MINPOS  IS  NOT  ALLOWED  AS  AN  A R 6 

00000  735 

tnz 

*•*■3 

IF  NOT  MINPOS*  GO  ON 

00000736 

LD  X7 

3, DU 

SFT  UNDERFLOW  FAULT 

OOOO  O 73  7 

TR  A 

8 ROUND 

GO  TO  SROUNDC.EU.)  WITH  JNFLO  JY  1 

0 0000  73  8 

FS  T 

ARG1 

STORE  IN  D0U3LE  PRECISION  AR E A 

COOOO 739 

ST  Z 

ARG1  -*-1 

CLEAR  FXTRA  PRECISION 

Cl  000 Q 74  0 

ea  xn 

3,1  • 

OOCjr  741 

FLD 

0/0 

00000742 

FCMP 

' MINN EG  ' 

~M  rNTlTG  mror  Allowed  AS  AN  ARGUMENT 

COBOL  74  5 

TNZ 

*♦2 

GO  ON  IF  ARG.  IS  NOT  MINN  EG 

C 0000  74  4 

TR  A 

BROUND 

AND  GO  TO  GROUND  PART  .60. 

0 COO 0 74  5 

FCMP 

MI NPOS 

MINPOS  IS  NOT  ALLOWED  AS  AN  ARGUMENT 

OOOOC  74  6 

TNZ 

* + 3 

GO  ON  IF  NOT  MINPOS 

00330074  7“ 

LD  XZ 

3,DU 

SET  UNDERFLOW  FAULT 

C 3 00  0 7 A 3 

' TR  A 

SPOON'D'  “ 

GO  TO  "P'FOUND  f . EJ'.r  w!Th  Eu  Sy  ONE 

00030  74  ? 

FST 

ARG2 

SAME  FOR  A R 62 

0 000  0 75  0 

STZ 

ARG7+T 

CLEAR  EXTRA  PRECISION. 

CCOUC 751 

I FE 

*1/2/4 

IF  BPACE8,  DO  NEXT  FOUR  LINES 

CUC00752 

LD  A 

0/G 

C 0030  753 

LDQ 

1 /O 

00000754 

STA 

A'R  6 T ' 

0(5c:r:  ’55 

STQ 

ARG1 ♦ 1 

OOOOU  756 

L D A 

*1  ,DL 

STORE  CALL  TYPE  IN  CALL.  CE9=2,  'others* 

1 COOOO  75  7 

STA 

CALL 

0 0000  75  3 

EN  DM 

ENTER 

00000759 

MACRO 

OOOOC  760 

STXO 

TEMP 

00000761 

LDXO 

COUNT 

GET  SHIFT  COUNT  AND  SKIP  IF  ZERO 

00000762 

r — 


i 

i 


1 

T Z F 

■ «u 

00000763 

CMPxO 

57,  OU 

SHIFT  NORMAL  IF  .LE.  36  BIT  SHIFT 

00000  764 

If. 

TM  1 

*♦11 

00000765 

CMPXO 

63, OU 

SHIFT  .GE.  63,  CHECK  FOR  ONES 

00000  766 

tmi 

* ♦ 6 

00000  767 

CANA 

= 077  77  777  777  77 

CHECK  FOR  ANY  ONES 

00000  766 

TZE 

* + 2 

If  NOT  THEN  SKIP,  ELSE  GO  ON 

00000769 

L 0 X 7 

1 ,0U 

SET  RESIDUE  INDICATOR 

00000770 

DFLD 

ZERO 

ALL  BITS  WERE  SHIFTED  OUT  SO  LOAD  ZERO 

OU000771 

TR  A 

*♦5 

NOW  PREPARE  TO  RETURN 

00000  772 

CANA 

RSTABL-37,0 

CHECK  FOR  ONES  THAT  WILL  BE  LOST 

00000773 

TZE 

* + 2 

IF  NONE  THEN  GO  SHIFT 

00000  774 

I'* 

L D X 7 

1 , DU 

ELSE  SET  RESIDUE  INDICATOR 

00000  775 

LRl 

0,0 

SHIFT  MANTISSA  RIGHT  COUNT  BITS 

00000  776 

LDXO 

TEMP 

00000777 

I 

ENOW 

ARSHF  T 

00000  77  3 

NEGU 

MACRO 

00000  779 

OF  ST 

SAFE 

00000  780 

DFLD 

R.U. 

COOOO  7 81 

f N EG 

00000782 

OF  ST 

R.U. 

OOOCO  783 

OFLO 

SA  FF 

00090784 

ENOW 

NEGU 

00000785 

8PACEB 

ENTER 

2 

ENTRY  TO  CONVERT  EXTENDED  TO  BPA 

C 0000  786 

OFLO 

ARG1 

00000737 

tr  a 

BROUND 

0000  0 788 

0P AMUL 

ENTER 

1 

00000789 

FLO 

ARG1 

MAKE  ARGS  POSITIVE  IF  NEC. 

00000  790 

tpl 

SPAM  . 1 

SET  SC  IF  WE  CHANGE  SIGN 

OCOOO  791 

FN  EG 

00000792 

FST 

ARG1 

00000793 

L0X6 

1 , DU 

00000  794 

8PAM. 1 

FLO 

A R G2 

00030795 

TPL 

BPA* . 3 

00000796 

. 

FNFG 

00000  797 

[ 

FST 

ARG2 

00000793 

CMPX6 

1 /Du 

NEED  TO  FLIP  SC. 

00000  799 

TZE 

BPAM . 2 

00000800 

LD  X6 

1 ,DU 

00000301 

TR  A 

8 P AM  . 3 

00000302 

BPA». 2 

LD  X6 

0 , 0 U 

0 003  0 8.3  3 

BPAW.3 

FL  0 

ARG1 

MULTIPLE,  DO  NOT  NORMALIZE 

00000804 

UFM 

ARG2 

OOOCO  305 

TEU 

PPAWEU 

00000806 

TEO 

QPAMtO 

OQOOC  807 

FNO 

EXPONENT  IS  OK/  TRY  TO  NORMALIZE. 

00000808 

H 

TEU 

BPAMEU 

WE  CAN  ONLY  WAKE  THE  EXPONENT  SMALLER. 

C 003 3 809 

CMPX  6 

0,  DU 

GET  SIGN  ANSWER  CORRECT. 

0000081 0 

I • 

TZE 

BRO'JND 

00000811 

p - - 

FN  EG 

0000081 2 

Tr  a 

BROUND 

OO0OO813 

BPAMEU 

FNO 

GET  THE  ANSWER  IN  THE  BEST  FORM  POSS. 

00000814 

1 


B5 


L X L4 

OPT! ON 

STZ 

BPAFLT 

CMPX  6" 

U,DU 

TZE 

.EU. 

FNEG 

TR  A 

• FU. 

a P A M E 0 

LXL4 

OPTION 

STZ 

BPAFLT 

DEST 

AN'S 

LDE 

0,  DU 

* NEGATE  BEFORE  eROUl 

CMPX6 

0 , D'J 

TZE 

BPAM.4 

FNEG 

BPAM.4 

FNC 

ADE 

ARG1 

ADE 

ARG2 

TEO 

.EO. 

TRA 

8 ROUND 

BPADIV 

ENTER 

1 

FLD  ~ 

argT 

tpl 

BPAV. 1 

FNEG 

D F ST 

A R G 1 

LDX6 

1 ,DU 

BPAV. 1 

FLD 

AR  G2 

TZE 

'ZEDDDV 

tpl 

BPAV. 3 

FNEG 

DF  ST 

ARG2 

CMPX  6 

T ,D'J 

TZE 

BPAV.  2 

L D X 6 

"TWOU 

TRA 

BPAV . 3 

BPAV. 2 

L0X6 

0,OU 

BPAV . 3 

STZ 

SHIFT 

FL  0 

ARG2 

sta 

TEMP 

FLD 

ARG1 

CMPA 

TEMP 

TIM  I 

BPAV. 5 

tnz 

*♦5 

CMPX  6 

C , DU 

TZE 

* + 3 

CMPX" 

=020001 

TZE 

‘ + 2 

ACS 

SHIFT 

LRS 

1 

TR  A 

BPAV  . 4 

BPAV. 5 

FLD 

ARG1 

BPAV. 4 

OVF 

_ TEMP"" 

CMPO 

0,0  J 

BEFORE  TRANSFER,  MAKE  SURE  ANS  IS  POS. 


oooo'osi  s 
0000081 6 
00000  81  7 
0000081  8 
OUUO'O  8 1 9 
00000  380 
'0  0000' "821 
0 000  0 8 2 2 


SHIFT  TO  NORMALIZE.  R(E)  WILL  HAVE  NE  G , GOOC3324 
IF  NEEDED  " 00C0C325 

00000826 

00000527 

0 COO  0 82  3 

SHIFT  NUM.  IT  MAY  BE  ENOUGH  TO  Fix  THE  CCOOO  S>9 

overflow,  this  first  add  can't  overflow.  00000830 
BECAUSE  BOTH  EXP.  OF  VAR.  ARE  POSITIVE  00000831 


SHIFT  NUM.  IT  MAY  BE  ENOUGH  TO  Fix  THE  CCOOO  3*9 

overflow,  this  first  add  can't  overflow.  00000830 
BECAUSE  BOTH'  EXP.  OF  VAR.  ARE  POSITIVE  00000831 

00000  832 

00000833 

00000834 

MAKE  BOTH  MANTISSA  POSITIVE  CCOOO  835 

0 0000  83  6 
OTjOTTC  837  " 
00000  838 

tnroinrsTv 

00000840 

— — ouoooBcr 

00000842 
CD000843 
30000844 
" 00000  845' 

00000346 

U'O'J'JJ  84  7' 

COOOO  348 

DD0D034  9 

COOOO  850 

'SEE  IF  DIVIDEND  rs  lA'RG'ER"  THAT  DTVT 53~R  OTJOOOT5T' 

THIS  WORKS  IF  BOTH  ARGS  ARE  NORMALIZED  00000852 


IF  THE  MANTISSAS  ARE  .5  AND  THE  SIGN 
WAS  CHANGED,  THEN  WE  MOST  NOT  ADD 
ONE  TO  THE  SHIFT  COUNT 


00000854 
~TrOOOTT555“ 
00000  856 
00000357 
00000353 
uujou  85  y 
00000860 
'00000361 
00000862 
CCOOO  86T 
OGOOC 864 
0 O'J'JO  36T 
COOOO  366 


F ILL  RESIDUE  INC  If 
SEE  IF  ANS  IS  NEG. 


CL£AJ?_  EXPQN£ilI_t  0__?R£  V ENT-  ERROR  IN  NEG 


0 .DU 

SHIFT 

0 . DU_  

EXP.  1 

O.DU  . _ . 

=0002000. DU 


=0002000. DU 

ANS 
A R G 1 

=0776000. DU 
28 


00000667 

00000869 
00000870 
0 000  J 87 1 
00000872 
00000873 
0 OOP  0 87  A 
00000875 
..O0QQ0  876 
0 000  J 87  7 

aOQflOSTE 

00000879 

00000879 

onnnnssn 
00000381 
— 00000  38  2- 
00000883 


LOG 

ARG2 

00000865 

ANQ 

=0776000. DU 

00000886 

GRS 

28 

00000387 

STQ 

TEMP 

00000385 

S d A 

TEMP 

00000889 

EDO 

ANS 

00000890 

ANQ 

=0776000. DU 

00000  891 

QRS 

2 3 . . 

00000892 

STQ 

TEMP 

UOOO J 393 

AOA 

TEMP 

OQOOD89A 

C M PA 

= 128 

00000395 

tpl 

EXPO 

CHECK  IF 

EXP 

IN  RANGE  00000896 

CMPA 

=-128 

OOCOO  897 

TMI 

EXPU 

CC00039S 

als 

28 

00000  899 

STA 

TEMP 

_ 00000900 

OFLD 

ANS 

00000901 

LDE 

TEMP 

00000902 

DF  ST 

ANS 

00000903 

CMPX  7 

1 ,ou 

I F THERE 

WAS 

A REMAINDER.  ADD  ONE  COOP 090 A 

TN2 

BR0UN9 

00000905 

STE 

APOONE 

00000  906 

FED 

AOOONE 

00000907 

TRA 

9R0UND 

00000908 

EXPO 

E X L A 

OPTION 

00000909 

S T 2 

BPAFLT 

CO.OOC  91 0 

TR  A 

• EO. 

0000091 1 

E XPU 

L X LA 

OPTION 

00000912 

ST2 

BPAFLT 

00000  91  3 

TRA 

.EU. 

00000  91  A 

2 ERODV 

EDA 

A.OL 

DIVISION 

9Y 

2ER0  CASE  00000915 

STA 

BPAFLT 

. ..  . ..  . 0000.0  916 

TRA 

EXIT 

0000091 7 

8P  A A 00 

ENTER 

1 

0000091 8 

FLO 
TM 
0 F ST 

SUBIN  Of  LD 


A RG2 



ARG1 

EXIT 

R.U. 

ARGt 

ADOI 

-8  .31. 

EXIT 
T ESI  1 
R.U. 

AO02 


pfld 

tpl 

FN  £ G 

OF  CMP 

TZ  E 

TM  I 

OFLD 

CFST 

OFLO 

TP_A 

OFLD 

OFCMP 

tpl 

-F.ti.EG  . 
NEGU 
L 0 X 6 
STE 

Loxn 


GET  ARG2,  TEST  FOR  ZERO,  STORE  IN  R.'J. 
ARG?  = 0,  LOAO  ARG1  AND  RETURN 


BPASUa  COMES  -MERE.  GET  ARGT. 

IF  A R G 1 * 0,  PREPARE  TO  LEAVE  BP  AA  00 

OUR  AN  SUER  IS  ARG2 

GO  TO  EXIT,  ONE  ARG  WAS  ZERO 

WE  WILL  MAKE  BOTH  POSITIVE  AND  COMPARE 

DfCMG  JQ-E-S  NQi_UQRK  CiiQU  . T9Z6J 

MAKE  2N0  ARG  POSITIVE  FOR  OFCMP 


MAKE  RCEAQ)  ♦ 

MAKE  THE  COMPARISON 

IF  EQUAL,  00  NOT  SWITCH 

- JLf.  RCEAQ)  .LT.  R.U-,  LEAVE 

RCEAQ)  MAY  8E  WRONG  SIGN,  SO  RELO 
AJJB_-SJV3ICil-JLC£AQ)_-AN-D  R.U._  


00000  VI  9 

—aooaoua- 

00000921 
00000  922 
OUOOU  92  3 
0000092 4 
OOOOG  92  5 

0000  3 926- 

00000927 
MD0QHS23 
00000929 
- 00003930 
00000931 
n nnn  ion 


00000933 
00000934 
00000  93  S 
_03303936 
00000937 


RECONSTITUTE  THE  SIGN  OF  RCEAQ) 
RCEAQ)  MUST  BE  SET  POSITIVE 


NEG.AI.E-  RCE AQ)  AND  R. 


-SXI-SXGN 


_1.,PU  _ 
SAVEXP 



Af.'XO 

=0  776000,  OU 

S6X0 

SAVEXP 

TOV 

*♦2 

TR  A 

* + 2 

LDXO 

=0200000, DU 

....  . TM! 

A 0 08  . 

TZE 

A008 

E A Q 

n,-0 

QRL 

10 

STQ 

COUNT 

LOQ 

= 0 

ARSHF1 

COUNT 

LDE 

R.U. 

DF  AO 

R.u. 

CMPX6 

0,DU 

TZE 

A001 1 . 

TEO 

*♦3 

. . LEU 

A0D9 

TR  A 

ADD1  0 

GET  EXPONENTS,  COMPUTE  SHIFT  COUN 


IF  OVERFLOW,  THE  SHIFT  NUMBER 

. W AS.  TOO  .BOG,  AND  SO  HE.  3.0  AD_  6 4 

AS  THE  SHIFT  COUNT 

IF  SHIFT-JS  NEGATIVE  NO  SHIFT  NEI 
NO  SHIFT  NEEDED,  SHIFT  COUNT  WAS 


STORE  SHIFT  COUNT  IN.C.OUNT 


LOAD  CORRECT  EXPONENT  FOR  RCEAQ) 

A DO  „TH  E TWO_.NU-MB.ER-S  AND  NORMALIZE. 
TEST  FOR  SIGN  CORRECTION 
AND  SKIP  TO  A 0.0 1 1 IF  NONE  NEEDED 
GO  TO  CORRECTION  FOR  EXPONENT  OFLC 

AND  AGAIN  IF  NO  FAULT  OCCURRED 
MINPOS  IS  THE  RESULT  0.F  ONE  OVER  F L 


AO  A R G 1 00000939 

00000940 

00000941 

00000947 

00000943 

QQQQQ  944- 

00003945 

C.HAN3E 000339.46. 

G JO J J 94  7 

130000343. 

IT  00000949 

Q000095Q- 

00000951 

30000952 

0 000  0 95  3 

-0  OOP  0 95  4 

00000955 

ZERO  00330957 
30000953 

00003959 

00000960 

00000961 

0000  0942 

0003  3 963 
. 00000  964 

00003965 
QQOQO  966 
0 00000967 

noooo  9&a 

00000969 
LOW  Q 000  0 97Q- 


J 


. 


* - 

CORRECT  NEGATIVE  IS  NINNEG#  BUT  NEGATIVE 
M1WPQS  GIVES  LIND  F R F i OU*  MTNMFT.  IS  NOT 

00000971 

ALLOWED  IN  INTERVAL/  SO  GIVE  IT  OVERFLOW 
TEU  TURN  OFF  UNOERFLOW  IF  WE  NEGATED 

00000973 

00000974 

FN£G 

TEU 

TURN 

MINPOS  . 

MLGATt  THE  MN.T1SSA 

00000975 
00000  976 

*♦1 

EOON 

TURN  EO  BACK  QN 

00000977 

0000097ft 

TRA 

A0D1  1 

WE  ARE  THROUGH  NEGATING  OVERFLOW 

M INN  EG  IS  THE  RESULT  OF  ONE  UNDERFLOW; 
CORRECT  NEGATIVE  IS  MINPOS/  BUT  NEGATIVE 
M INN  E G GIVES  OVERFtOWj  MLNPQJ  15  NQT 

00000979 
OO0O0  980 
00000  9S1 
£0000942 

ALLOWED  IN  INTERVAL/  SO  GIVE  IT  UNFlOW. 
TEO  TURNS  THE  OVERFLOW  INDICATOR  OFF  IF 

00000983 
00000  984 

WE  NEGATED  MINNEG. 

00300985 

ADD9 

FNEG 

NEGATE  the  UNDERFLOW  ANSWER  . 

_0 0000  946 

TEO 

*♦1 

TURN  OFLO  OFF  IF  IT  IS  ON 

00000987 

TURN 

EUON 

_T  UR N_  EU  INDICATOR  BACK  ON 

00000  988 

TRA 

*001  1 

AND  GO  TO  ADD11 

C 0000  989 

*0010 

FNFG 

NO  PROBLEM,  j UST  NES A TE 

00300  993 

A D 0 1 1 

CWPX7 

0#  DU 

CHECK  RESIDUE  INDICATOR/  EXIT  IF  ZERO 

0C900991 

tnz 

A 0 0 1 ? 

00000  79? 

DF  CMP 

ZERO 

EXIT  IF  ZERO.  ELSE  GO  TO  8R0UND 

00000993 

TZE 

EXIT 

00900  994 

TRA 

8R0UND 

OCOOC  995 

00333996 

00300997 

AT  THIS 

POINT  WE 

MAY  NEED  A RESIDUE  CORRECTION  . 

00300994. 

WE  HAVE 

THE  FOUR 

FOLLOWING  POSSIBLE  CASES.. 

COOOO  999 

1,  sc  = 

0 AND  RCEAQ)  .GE.  0 

C0901  300. 

?.  sc  = 

0 AND  RCEAQ)  .LT.  0 

0000 1 001 

3.  SC  = 

1 AND  RCEAQ)  .GE.  0 

00001 30? 

4.  SC  = 

1 AND  RCEAQ)  .LT.  0 

00001335 

00001 304 

CASES 

i 

AND  ?. 

IF  THE  EXTRA  PRECISION  WORD 

0CC01  OC  5 

IS  0/  SET  IT  TO  1_. 

00001  00.6 

AND  THEN  GO  TO  BROUN  D . 

00031 007 

CASES 

T 

AND  4. 

IF  THE  EXTRA  PRECISION  WORD 

00001 303 

IS  0/  SUBTRACT  1 FROM  THE 

CC001 309 

FULL  MANTISSA/  AND  THEN  GO 

00001 010 

TO  BROUNO. 

00001  01 1 

0000’  01  ? 

03001013 

*001? 

CMP*  6 

0/  ou 

IDENTIFY  THE  CASE.  CHECK  SC 

00001  314 

TNZ 

A 001  3 

IF  SC  = 1 / GO  TO  A0D1 3 

00001 31 5 

DFST 

TEMP 

00001 016 

LOQ 

TEfcP+1 

00001 01 7 

TNZ 

B ROUND 

IF  Q NOT  ZERO/  WE  ARE  DONE  HERE 

00001  01  s 

LOQ 

= 1 

00301  019 

S TO 

TEMP+1 

SET  Q TO  1 REPLACE  IN  TEMP 

30001 J?0 

DFLO 

TEMP 

LOAD  TEMP  AND  GO  TO  BROUND 

00001 021 

TRA 

3R01JND 

00001 022 

I 

! 

I 


■i 


1 


*0013  OFST 

T€MP 

WE  HAVE  CASE  3 OR  4 

00001 023 

. 1 QQ 

TF.MP 

rHFf*  fns  an n/fbo  aaia  i if  *n 

nflOQ i Q pl 

TNZ 

BROUND 

00001 025 

OFLO 

TEMP 

ELSE  SUBTRACI  ONE  FROM  AQ 

00001026 

STE 

SUdONE 

00001 022 

if  S3 

SUdONE  . . .. 

tiudoioza 

TR  A 

BROUND 

0000102V 

-fiPRSUR  . E N TRR 

...  .1  . 

nonn i n 

DFLD 

ARG2 

LOAD  ARG2,  TEST  FOR  ZERO.  NEGATE 

00001031 

TNZ 

*♦3  

-OQJO  l U3 Z 

FLO 

ARG1 

00001 033 

TRA 

EXIT  

000C1034 

FNEG 

NEGATE  ANO  GO  TO  SU3IN  IN  BP  A A 0 0 

00001 03S 

OFST 

R.U. 

DMOI'Hl, 

OFST 

A R b2 

00001037 

TRA 

SUBIN 

......  . . . _ . _ 

00001 038 

8R0UND  FNO 

NORMALIZE  THE  R(EAQ)  IF  IT  IS  NOT 

00001 039 

. S T.Z  _ 

BPAFLT 

PREPAfLE  PEI  URN.  \LAK1  API  F _ 

nnnn  i run 

OF  CMP 

ZERO 

IF  ANSWER  IS  ZERO  WE  ARE  DONE 

00001 041 

TZf 

. . i .XJ  T 

THEN  GO  TO  EXIT, WE  ARE  DONE 

00301  04z? 

LXL4 

OPTI ON 

00001043 

. TUI 

* *4 

. ODQOJ-04-4- . 

TZE 

* +3 

000C1 045 

CMPX4 

6-,BU- 

_ 0000.1-046- 

TM  I 

*♦2 

000C1 047 

LDX4 

4.  DU 

00001 046 

TEO 

• EO. 

DID  OP  CAUSE  EXP  OVERFLOW  OR  UNDERF. 

00001 04V 

TE'J 

» EIJ  * 

__  nnpoi  nsp 

OFCMP 

MI  NPOS 

IF  ANSWER  IS  MINPOS,  GO  TO  .EU.1 

00001 051 

tze 

. E 1 

nnnn ins? 

OFCMP 

NEGBND 

IF  ANS . >*NE GBND^NO  ERROR 

00001053 

T M I 

.EO. 

00001 054 

OFCMP 

BIGC  HK 

IF  BIGGER, THEN  OELOW  EVEN  THOUGH 

00001 055 

TM1 

*♦3 

E.T_IS  .SMALLER  THAN  MAXPOS  _ _ . 

00001056 

TZE 

BRNO . 1 

GO  LOAD  OKTAB  JUMP  VALUE  . 

00001 C57 

...  TRA 

«.£  d,  .. 

UF  HAWF  QVPRFLQU#  U.1  T H Q Lil  EQ_  ON.  . 

0000-1-0511- 

C A N A Q 

CHKDBL 

IF  NO  ONES  IN  08 L PREC  PART  OF  MANT 

00001 059 

T^E_. 

.6  XJ.T 

THEN  WE  NEED  NO  BROUNDING 

00001 060 

OFCMP 

ZERO 

CHOOSE  UPPER  OR  LOWER  HALF  OF  WORD 

00001 061 

TM1 

* + 3 

T_AK_£  LOWFR  HA1  F IF  MI  NUS 

. _QQ001  062 

BRND.1  LDX4 

OK  T A B * 4 

GET  JUMP  VALUE  FROM  TABLE 

00001 063 

TRA 

NO  F l T * 4 

GO  TO  RIGHT  3R0UNDING 

000.0-1  064. 

L X L4 

0KTAB,4 

GET  JUMP  VALUE  FROM  TABLE 

0000106? 

TRA 

NO  FL T # 4 

GO  TO  BROUNDING  FOR  MpigS 

00001 C66 

.EO.  FCMP 

ZERO 

TEST  WHETHER  POSITIVE  OR  NEGATIVE 

00001 067 

TM  I 

.EO.  1 

IF  NEGATIVE  GO  TO  .EQ.1 

0 0 00.10  6 8 

L D X 2 

E 0 T A B#  4 

LOAD  FAULT  VECTOR  FROM  TABLE 

0003 1 069 

FLO 

MAXPOS 

LOAD  POSITIVE  ANSWER 

00001  07Q. 

TRA 

EXIT 

00001 071 

• EO. 1 FLO 

NEGBND 

LOAD  NFGATIVE  ANSWER 

00001 072 

LXL2 

E OT  * 8/ 4 

LOAD  FAULT  VECTOR 

00001 073 

TRA 

EXIT 

WE  ARE.  THROUGH - - - 

QQQ0JO74- 

BIO 





tU  LDX2 

5 /DU 

LOAD  UNDERFLOW 

FOR  ALL  OPTIONS 

00001 075 

FCMP 

2 £80 

_ I ESI  if  PQS_III 

V£  OR  NE  .'.A  I J u£ 

nnnn  i07f> 

TM  I 

.EU.  2 

I F NEGATIVE/  GO  TO  .EU.2 

00001 077 

CMPX  U 

A/DU 

ROUNDING  TAXES 

SPECIAL  TREATMENT 

00001078 

tie 

*♦3 

SO  SKIP  THESE 

TWO  LINES 

C0001 079 

ELD 

ELI  A 0E>  / 4 

LOAD  CORRECT  answer  FOR  U/L/T/A 

- . uoaaiwio 

TRA 

EXIT 

and  GO  TO  EXIT 

C0J01 081 

STE 

SAVEXP 

C HEC  K IF  UNFLO 

BY  1 

0 0001  082 

IDXO 

SAVEXP 

R(EAQ)  NOW  HAS 

ZERO  FOR  R OPTION 

C0001 085 

CMPXO 

EU  ON  E 

IF  UNFLO  BY  1 , 

LOAD  posbnd 

00001034 

TNZ 

♦ ♦3 

00001 J85 

. Ell.  1 FID 

POSRND 

00001 086 

TRA 

EXIT 

WE  ARE  DONE  WITH  POSITIVE  CASE 

00001 087 

ELD 

..  JJsJLCL  _ 

IF  EU  BY  >1  LOAD  ZERO 

ODQD 1 0 HA 

TRA  EXIT  AND  GO  TO  EXIT  00001089 


.Ell. 2 CMPX.4  _ LE  NOT  RQU-ALD.  _LDAD  V Ai.0  E f JLQM_  I ILbLE ^QOOQIMO- 


T2E 

FLO 

TRA 

STE 

* ♦ 3 

EllTABN./A 

LOAD  RIGHT  ANSWER  FOR  U/L/I/A 

0 GOO  1 09 1 
_ HQQai-092- 

EXIT 

SAVEXP 

AND  go  to  EXIT 

CHECK  If  UNFLO  0Y  1 

00001093 
00001 094 

LDXO 

CMPXO 

T N 2 

ELD 

SAVEXP 

EUONE 

R(EAQ)  NOW  HAS  ZERO  FOR  R OPTION 

IF  UNFLO  BY  t/  LOAD  MAXNEG 

00001 095 
_OIiQOJ  096 

* + 3 

MAXNEG 

ELSE  GO  LOAD  ZERO 

00001097 
. Q0MU93 

*- 

TRA 

fid. 

TRA 

EXIT 

ZERO 

LOAD  ZERO  AS  THE  ANSWER/  EU  > 1 

00001 099 

EXIT 

AND  GO  TO  EXIT 

JHE  URSX  BRQUNDI  NG_  k£  DO  IS  U(+,~), 

00001101 
00001  102 



T <- ) , AND  AC*).  IN  ALL  CASES/  WE 

ADO  UTA,  WHICH  GENERATES  A CARRY 

INTO  THE  SINGLE  PREC.  ANSWER  IF  THERE 
ARE  ANY  1*S  IN  THE  EXTRA  PREC. WORD 

C0001 103 

P.Q0Q1  m 

00001 105 
00001 106 

noflt 

SEE 

uta 

00001 107 

DF  AD 

L'TA 

JIQQQJ  ms  _ 

TEO 

* + 3 

OFLO  IS  POSITIVE.  HAVE  U OR  A 

00001109 

TEU 

*♦■5 

UNDERFLOW  IS  NEGATIVE^,  T 

.aaQajjj.a_ 

TRA 

EXIT 

GO  TO  EXIT  IF  ANSWER  OX 

00001  111 

LD  X2 

2/DU 

SET  INFINITY  FOR  OVERFLOW 

00001  112 

DF  LD 

MAXPOS 

WANT  MAXPOS  AS  ANSWE? 

00001 113 

TRA 

EXIT 

00001  114 

L D X 2 

3/OU 

SET  UNDERFLOW 

00001  115 

DFLD 

7E  RO 

LOAD  ZERO  ON  T OPTION 

COO01  116 

TRA 

EXIT 

E ND  OF  UTA  BROUN  DING 

00001  11  7 

* 

ROUND  POSITIVE  NUMBER.  BIT  28  MUST 

00001  11  8 

* 

BE  1 TO  GO  UP/  ELSE  CHOP  OFF  EXTRA 

00001  119 

• 

PRECISION. 

C0JJ01  120 

RPOS 

STE 

POSRND 

U0001 121 

DEAD 

POSRND 

00001 122 

TEO 

* + 2 

OVERFLOW  IS  POSITIVE 

00001 123 

TRA 

EXIT 

00001124 

LDX2 

2 /DU 

SET  INFINITY 

00001 125 

DFLD 

MAXPOS 

WANT  MAXPOS  AS  ANSWER 

00001 126 

Bll 


TR  A 

F X I T 

DfNJA/p  A *J££ATjyf  AZJftffO  ftQUNO 

00001 127 
C0Q01  13*- 

* 

* 

value  in  middle  toward  zero. 

BIT  2b  MUST  BE-  0 TO  GQ  AWAY 

0002 1129 
.000111  liO 

* 

RNE  G 

S T E 

NEG8ND 

FROM  ZERO. 

00001  151 
0000 1 132 

OF  AD 

T.LU- 

NEGRNO 

*±z 

UNDFRFLQU  IS  NF  G A TIM? 

00001133 
nnnm  i 3 u 

TR  A 

LDJ(2 

EXIT 

3/  QJ 

S£L_UZZD£RFL04^1iItFLO-  AT  1 

00031135 
OOQQU-SA 

EXIT 

OFLD 

SXL2 

MAXNEG 

9PAFLT 

SO  LOAD  MAXNEG/  NOT  0 

&EIUSN  / AULT  VEC  IOfi.  IN  COMMON-  . 

00001 137 
nnnm  i is 

LXL? 

ryox? 

CALL 

00001 139 
nnnm  iz.n 

T Z t 

EAXO 

FST 

TR  A 

EXIT.1 

4/1* 

O/O 

EXIT. 2 

AND  BROUNDED  ANSWER  IN  ARGUMENT 

00001 141 
0003  1142 
00001 143 
OQOQ1 1 44 

EXIT.1 

f A X r' 

FST 

3*1* 

0*0 

00001 145 
nnnm  14a 

EXIT. 2 

LREG 

Rtl 

. RSS  A 

0003 1 147 

nnnm  na. 

END 

00031 149 
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